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A  COMPUTER  PROGRAM  FOR  OPTIMAL  CONTROL  OF  WATER  SUPPLY 
PUMP  STATIONS:  DEVELOPMENT  AND  TESTING 


1  INTRODUCTION 


Background 

Energy  conservation  in  the  United  States  has  been  an  important  issue  ever  since  the  oil  embargo  of 
1973  when  the  public  became  very  much  aware  of  the  unstable  and  limited  availability  of  energy- 
producing  resources.  Shortly  after  the  oil  embargo,  the  attention  of  the  news  media,  political  leaders,  and 
the  American  public  was  focused  on  energy-related  topics  which,  at  the  time,  was  rare-especially  for 
technical  issues.1  Even  though  the  procurement  of  energy-producing  resources  has  stabilized  somewhat, 
supplies  are  still  limited  and  the  cost  of  these  resources  continues  to  rise. 

Water  utilities  are  major  consumers  of  electrical  energy.  In  fact,  up  to  7  percent  of  the  electricity 
consumed  annually  in  the  United  States  is  used  by  the  municipal  water  works  industry.2  Pumping  of 
water  comprises  the  major  fraction  of  electrical  usage  at  water  utilities,  accounting  for  as  much  as  80  to 
95  percent  of  the  energy  requirement  in  some  systems.3  Furthennorc,  energy  costs  are  the  most  expensive 
single  item  in  the  operating  budgets  of  many  water  utilities.4 

The  same  situation  exists  at  many  Army  installations  that  have  their  own  water  supply/distribution 
systems.  These  energy  costs  are  especially  devastating  in  light  of  the  continual  budget  cuts  within  the 
military  services.  Moreover,  the  Army  has  directed  its  elements  to  explore  every  feasible  option  for 
reducing  energy  consumption.5 

It  is  apparent  that  a  major  effort  to  reduce  the  electrical  consumption  and  associated  cost  at  water 
utilities  should  focus  on  efficient  operation  of  raw  and  finished  water  pumps.  This  reduction  could  be 
accomplished  by  optimal  real-time  control  of  water  supply  pump  stations.  Optimal  control  strategies  for 
water  distribution  system  operation  could  be  formulated  using  powerful  optimization  techniques  (Chapter 
2  presents  a  literature  survey  on  this  technology.)  To  be  effective,  such  a  method  would  allow  operators 
to  schedule  pump  operations  such  that  the  electrical  consumption  of  the  entire  pump  station  is  minimized 
while  adequate  system  storage  for  fire  protection  and  sufficient  system  pressures  are  maintained. 

Optimal  control  of  water  supply  systems  may  be  welcome  to  plant  operators  in  both  the  Government 
and  private  sector.  According  to  a  1984  survey  by  an  American  Water  Works  Association  (AWWA) 
committee,  many  utilities  would  like  to  see  a  progression  towards  computer-assisted  or  even  computer- 
controlled  operation  of  water  systems.6  The  desire  for  such  a  movement  is  due  in  part  to  (1)  the 


1  D.  A.  Dreyfus,  "Persistent  National  Energy  Issue,"  ASCE  Journal  of  Professional  Issues  in  Engineering,  Vol  114,  No.  1  (1988), 

pp  18. 

2  D.  Brailey  and  A.  Jacobs,  "Energy  Management  in  the  Water  Works  Industry,"  Journal  of  the  New  England  Water  Works 
Association.  Vol  94.  No.  3  (.1980),  pp  216-239. 

5  J.  L.  Patton  and  M.  B.  Horsley,  "Curbing  the  Distribution  Energy  Appetite,"  Journal  of  the  American  Water  Works  Association 
(JAWWA),  Vol  77.  No.  6  (1980),  pp  314-320;  H.  F.  Kellis  and  M.  K.  Griffin,  "Energy  Costs  Reduction  Through  Operational 
Practices,"  Proceedings  of  the  AWWA  Annual  Convention,  Dallas,  TX  (American  Water  Works  Association  [AWWA],  1984. 

4  w  H.  Clingenpeel,  "Optimizing  Pump  Operating  Costs,"  JAWWA,  Vol  75,  No.  10  (1983),  pp  502-509;  T.  Clark,  ''Reducing 
Power  Costs  for  Pumping  Water,"  Opflow,  Vol  13,  No.  10  (1987). 

'  Army  Regulation  (AR)  11-27,  The  Army  Energy  Plan  (July  1989). 

*  AWWA  Computer  Assisted  Design  of  Water  Systems  Committee,  "Network  Analysis  Survey  1984,"  Proceedings  of  the  AWWA 
Annual  Conference,  Dallas,  TX  (1984). 
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increased  complexity  of  distribution  systems,  (2)  high  operating  costs,  (3)  computer  hardware  which  is 
becoming  cheaper,  faster,  and  more  readily  available,  (4)  managers  and  operators  who  are  less  threatened 
by  computers  than  in  past  years,  and  (5)  more  sophisticated  simulation  models,  optimization  algorithms, 
and  control  software.7 


Objective 

The  objective  of  this  work  was  to  develop  and  test  a  computerized  method  that  will  enable  Army 
water  utility  managers  and  operators  to  select  an  optimal  control  strategy  for  water  supply  pump  stations. 


Approach 

To  develop  the  optimal  control  system,  the  literature  was  surveyed  for  existing  technologies  that 
might  be  useful  (Chapter  2).  Two  that  appeared  promising  (a  nonlinear  optimization  algorithm  and  a 
simulation  routine)  were  combined  to  form  the  model  (Appendix  A).  Two  subroutines  were  written  next: 
one  for  data  entry  and  another  to  provide  computational  features  (Appendices  B  and  C,  respectively).  A 
second  subroutine  for  data  entry  was  also  written.  The  resulting  program  is  called  Optimal  Control  of 
Pumping  Stations  (OCOPS).  It  is  written  in  FORTRAN  77  and  can  be  used  on  any  IBM-compatible 
personal  computer  (PC). 

To  determine  this  strategy’s  effectiveness,  OCOPS  was  tested  at  the  Fort  Hood,  TX  water 
distribution  system.  The  results  were  analyzed  in  terms  of  OCOPS’  ability  to  reduce  energy  costs  by 
suggesting  schedules  for  pumping  at  nonpeak  periods. 


Scope 

The  program  developed  can  be  applied  to  any  distribution  system  and  is  not  limited  to  any  number 
of  pumps,  pump  stations,  or  tanks  within  the  system.  It  provides  an  optimal  control  strategy  based  on 
electricity  unit  charges,  but  can  be  modified  to  account  for  electrical  demand  charges  as  well  by  including 
additional  constraints. 


Mode  of  Technology  Transfer 

It  is  recommended  that  the  information  in  this  report  be  transferred  using  both  an  information 
exchange  and  implementation  package.  The  information  exchange  would  consist  of: 

1.  Briefings  on  the  program  and  reports  in  field  experience  to  be  reported  at  the  PROSPECT  course 
entitled  Water  Supply  Design  and  Rehabilitation,  which  is  widely  attended  by  DEH,  MACOM,  and 
Corps  District  personnel,  and 

2.  The  DEH  Digest. 


U.  Shamir,  "Compuicr  Applications  for  Real-Time  Operation  of  Water  Distribution  Systems."  Proceedings  of  the  ASCE  Specialty 
Conference  on  Computer  Applications  in  Water  Resources,  H.  C.  Tomo  (Ed.)  (American  Society  of  Civil  Engineers  [ASCE], 


19K5). 
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The  implementation  package  would  include: 


1.  A  fact  sheet  lor  distribution  at  the  DEI1  conference  and  a  recommendation  that  interested  DEH 
utility  personnel  attend  the  PROSPECT  course  entitled  Water  Supply  Design  and  Rehabilitation. 

2.  A  computer  program  user  guide  including  specifications  for  gathering  data  for  program 
inplcmcntation. 

It  is  also  recommended  that  this  computer  program  be  further  demonstrated  through  either  the  Facility 
Engineering  Application  Program  (FEAP),  or  via  a  Technology  Transfer  Test  Bed  (T3B). 

The  proponent  agency  for  this  computer  program  is  the  U.S.  Army  Engineering  and  Housing  Support 
Center  (USAEHSC).  The  POC  at  USAEHSC  is  the  Sanitary  Branch,  at  (703)  355-7963.  Technical  details 
on  the  use  of  this  computer  program  can  also  be  obtained  from  the  U.S.  Army  Engineer  Waterways 
Experiment  Station  (USAWES),  phone  (601)  634-2879  or  toll  free  at  (800)  522-6937,  extension  2879. 
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2  SYSTEM  DEVELOPMENT 


Literature  Review 

Much  of  the  previous  research  on  reducing  electrical  consumption  or  increasing  the  efficiency  of 
pump  operation  in  water  distribution  systems  has  relied  on  heuristic  approaches.  Such  general  measures 
to  improve  energy  efficiency  include  the  use  of  variable-speed  pumps,  pumping  during  periods  of  low 
electrical  costs,  use  of  diesel-powered  pumps  during  peak  demand  periods,  improving  driver  power  factors 
by  installing  capacitors,  reducing  the  number  of  pump  starts  and  shutdowns,  and  reducing  pump  heads  by 
cleaning  and  lining  mains  and  opening  partially  closed  valves.8 

Other  researchers  have  employed  mathematical  optimization  techniques  in  an  effort  to  develop 
optimal  control  strategies.  Most  of  this  effort  has  focused  on  using  dynamic  programming  since  storage 
tank  operation  is  well  suited  to  state-stage  discretization.9 

The  efficiency  of  dynamic  programming  deteriorates  rapidly  with  an  increase  in  the  number  of 
storage  tanks  in  the  distribution  system.  Some  researchers  have  suggested  decomposing  the  system  into 
subsystems  with  one  or  two  storage  tanks  and  solving  the  optimal  control  program  using  dynamic  pro¬ 
gramming.10 

Other  investigators  have  used  nonlinear  optimization  techniques  to  solve  the  optimal  control  problem. 
These  techniques  are  applied  to  simplified  representations  of  the  water  distribution  system’s  dynamics  and 
are  useful  with  highly  dimensional  systems.'1  Cohen  provides  a  very  good  summary  of  the  methods  and 
algorithms  developed  for  optimal  control  of  water  supply  pump  stations.12 

Very  little  work  has  been  done  to  integrate  optimization  algorithms  with  water  distribution  analysis 


*  G.  A.  Aldworth,  "Energy-Saving  Pump  Selection,"  JAWWA,  Vol  75,  No.  10  (1983),  pp  496-501;  D.  Brailey  and  A.  Jacobs;  J. 
Chao,  "Can  Off  Peak  Pumping  Cut  Utility  Power  Costs,"  JAWWA,  Vol  71,  No.  5  (1979),  pp  259-263;  T.  Clark;  W.  H. 
Clingenpcel;  G.  W.  Lackowitz  and  P.  J.  Pretrctti,  "Improving  Energy  Efficiency  Through  Computer  Modeling,"  JAWWA,  Vol 
75,  No.  10  (1983),  pp  510-515;  J.  L.  Patton  and  M.  B.  Horsley;  H.  F.  Rehcis  and  M.  K.  Griffin;  B.  S.  Aptowicz  et  al.,  "Using 
Elevated  Storage  and  Off-Peak  Pumping  to  Control  Energy  Costs,"  JAWWA,  Vol  79,  No.  1 1  (1987). 

9  B.  Coulbeck  and  C.  H.  Orr,  "Optimized  Pumping  in  Water  Supply  Systems,"  Proceedings  of  IFAC,  IX  Triennial  World 
Congress,  Budapest.  Hungary,  Vol  6  (1984);  L.  W.  Ormsbec,  et  al..  Techniques  for  Improving  Energy  Efficiency  at  Water 
Supply  Pumping  Stations,  Technical  Report  EL-87-16  (U.S.  Army  Waterways  Experiment  Station,  1987);  M.  H.  Sabel  and  O. 
J.  Helweg,  "Cost  Effective  Operation  of  Urban  Water  Supply  System  Using  Dynamic  Programming,"  Wc  er  Resources  Bulletin, 
Vol  21,  No.  1  (1985),  pp  75-81;  U.  Shamir;  J.  L.  Solanas  and  M.  Verges,  "Extension  of  the  Dynamic  Programming  Successive 
Approximation  and  Its  Application  to  Automatic  Operational  Control  of  Water  Distribution  Systems,"  IFACllFORS  Symposium, 
Varna,  Bulgaria  (1974);  M.  J.  H.  Sterling  and  B.  Coulbeck,  "A  Dynamic  Programming  Solution  of  Optimization  of  Pumping 
Costs."  Proceedings  of  the  Institute  of  Civil  Engineers,  Vol  59  (1975),  pp  813-818;  A.  J.  Tarquin  and  J.  Dowdy,  "Optimal  Pump 
Operation  in  Water  Distribution,"  ASCE  Journal  of  Hydraulic  Engineering,  Vol  1 15,  No.  2  (1989),  pp  158-168. 

in  B.  Coulbeck  and  C.  H.  Orr,  "Optimized  Pumping  in  Water  Supply  Systems,"  Proceedings  of  the  3rd  JFAC  Symposium  on 
Control  of  Distributed  Parameter  Systems,  Toulouse,  France  (1982). 

"  Ccmbrano,  et  al.,  "Optimization  of  a  Multi-Reservoir  Water  Network  Using  a  Conjugate  Gradient  Technique.  A  Case  Study." 
8th  International  Conference  on  Analysis  and  Optimization  of  Systems  (Inslitut  National  de  Recherche  en  Informalique  et  en 
Auiomatique,  Antibes,  France,  1988);  B.  Coulbeck,  "Optimal  Operations  in  Non-Linear  Water  Networks,"  Optimal  Control 
Applications  and  Methods,  Vol  1  (1980),  pp  131-141;  B.  Coulbeck  and  M.  J.  H.  Sterling,  "Optimized  Control  of  Water 
Distribution  Systems,"  Proceedings  of  IEE,  Vol  125,  No.  9  (Institute  of  Electrical  Engineers  [IEE],  1978),  pp  1039-1044;  F. 
Fallsidc  and  P.  F.  Perry,  "Hierarchical  Optimization  of  a  Water-Supply  Network,"  Proceedings  of  IEE,  Vol  122,  No.  2  (1975), 
pp  202-208;  M.  J.  H.  Sterling  and  B.  Coulbeck,  "Optimization  of  Water  Pumping  Costs  by  Hierarchical  Methods,"  Proceedings 
of  the  Institute  of  Civil  Engineers,  Vol  59  (1975),  pp  789-797. 

12  G.  Cohen,  "Optimal  Control  of  Water  Supply  Networks,  Optimization  and  Control  of  Dynamic  Operational  Research  Models, 
S.  G.  Tzafetas  (Ed.)  (North  Holland  Publishing  Co.,  1982),  pp  251-276. 
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packages.  Although  a  few  models  have  beer  developed  using  this  procedure,13  they  consider  only 
optimal  design  and  not  optimal  operation. 


Developmental  Approach 

The  optimal  control  method  developed  in  this  study  is  conceptually  very  simple.  A  nonlinear 
optimization  algorithm,  GRG2,14  was  coupled  with  a  water  distribution  simulation  routine,  WADISO.15 
The  decision  variable,  or  quantity  that  can  be  changed  to  achieve  least-cos^  pumping,  is  the  fraction  of 
time  a  pump  is  operated  during  a  given  time  interval.  The  nonlinear  optimizer  and  the  simulation  routine 
were  linked  together  by  a  subroutine  called  GCOMP  which  computes  the  cost  of  operation  and  updates 
all  tank  water  levels. 

GRG2  systematically  adjusts  the  pump  operating  period  during  a  given  time  interval  and  provides 
this  information  to  the  GCOMP  routine.  Based  on  this  information,  GCOMP  finds  the  boundary 
conditions  and  passes  this  information  to  WADISO,  which  solves  the  equations  describing  the  pressure 
and  flow  distribution  in  a  water  system.  WADISO  returns  values  for  pump  head  and  discharge  for  each 
operating  pump,  flow  in  lines  connecting  tanks  to  the  system,  and  pressure  at  critical  nodes  to  the  GCOMP 
subroutine.  GCOMP  computes  the  electrical  consumption  and  resulting  operating  cost  for  the  specified 
time  period;  it  also  updates  tank  levels.  This  information  is  returned  to  GRG2  where  it  is  used  to  generate 
an  improved  operating  strategy,  and  the  process  is  repeated  until  an  optimal  solution  is  reached.  Figure 
1  illustrates  this  process. 

To  test  OCOPS,  an  experimental  run  was  to  be  conducted  at  Fort  Hood,  TX.  Chapters  3  and  4 
describe  the  Fort  Hood  water  distribution  system  and  provide  details  on  the  computer  model  and  the  test 
application.  Appendix  A  contains  a  detailed  description  of  the  solution  methodology;  GRG2,  WADISO, 
and  OCOPS  also  are  explained  further  in  this  section. 


"  L.  Onnsbee  and  D.  Contractor,  "Optimization  of  Hydraulic  Networks,"  International  Symposium  on  Urban  Hydrology, 
Hydraulics,  and  Sediment  Control,  Lexington,  KY  (1981 );  K.  E.  Lansey,  Optimal  Design  of  Large  Scale  Water  Distribution 
Systems  V  ter  Multiple  Loading  Conditions,  Ph.D.  Dissertation  (University  of  Texas  at  Austin,  1987), 

"  l,.  S.  Lasdon  and  A.  D.  Warcn,  GRG2  User's  Guide  (University  of  Texas,  1986). 

”  J.  Gcsslcr  and  T.  M.  Walski,  Water  Distribution  System  Optimization,  Technical  Report  EL-85-11  (U.S.  Army  Engineer 
Waterways  Experiment  Station.  1985). 
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Figure  1.  Optimal  controt  process. 
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3  DESCRIPTION  OF  THE  FORT  HOOD  WATER  DISTRIBUTION  SYSTEM 


Fori  Hood  is  located  in  cast  central  Texas  approximately  60  miles  north  of  Austin.  The  installation 
serves  a  population  of  more  than  64,000  and  is  one  of  the  largest  military  installations  in  the  United 
States.  The  post  is  divided  into  three  general  areas:  North  Fort  Hood,  training  areas,  and  the  Main  Post. 
The  Main  Post  contains  the  major  fraction  of  total  population,  building  space,  and  installation  utilities. 
Therefore,  only  the  distribution  system  serving  the  Main  Post  was  included  in  the  analysis. 


Water  Source 

Water  is  purchased  from  the  Bell  County  Water  Control  and  Improvement  District  (WCID)  and  is 
delivered  to  Fort  Hood  through  a  storage  and  interconnection  facility  adjacent  to  the  installation’s  main 
pump  station.  Two  5  million  gallon  (MG)*  ground  storage  tanks  at  the  facility  supply  water  to  Fort 
Hood.  Water  flows  from  the  storage  tanks  through  metered  interconnections  to  the  main  pump  station  and 
also  to  three  ground  storage  tanks  located  next  to  the  pump  station. 


Main  Pump  Station 

Water  is  delivered  to  the  Main  Post  through  five  vertical  turbine  pumps.  Table  1  lists  the  rated 
characteristics  of  each  pump.  Pump  operation  is  controlled  automatically  by  the  water  level  in  tank  5  in 
the  western  part  of  the  main  cantonment  area.  Sec  Table  2  for  level  settings  of  each  pump. 


Tabic  1 

Fort  Hood  Main  Pump  Station  Pump  Characteristics 


Pump  No. 

Serial  No. 

Type 

Rated 

Horsepower 

(hp) 

Rated 

Head 

(ft) 

Rated 

Discharge 

(gpm) 

1 

Verti-Line  D06876 

Vertical  Turbine 

100 

262 

1,150 

-*> 

4. 

Vcrti-Line  D06760 

Vertical  Turbine 

200 

262 

2,350 

3 

Fairbanks- Morse 

P2E3238 

Vertical  Turbine 

300 

262 

3,500 

4 

771H0437  270 

Vertical  Turbine 

500 

262 

6,000 

5 

Verti-Line  101449 

Vertical  Turbine 

500 

262 

6,000 

‘A  metric  conversion  table  appears  on  p  43. 
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The  values  shown  in  Table  2  represent  the  summer  pumping  policy.  This  pumping  sequence  results 
in  the  200-hp  pump  starting  first,  followed  by  the  300-hp,  100-hp,  new  500-hp  (pump  5),  and  old  500-hp 
(pump  4),  respectively.  The  winter  pumping  sequence  is  the  same  as  the  summer’s;  however,  the  tank 
elevations  arc  slightly  different.  Basically,  the  pumps  are  operated  so  as  to  keep  all  storage  tanks  as  full 
as  possible.  Figure  2  is  a  detail  of  the  main  pump  station  and  WCID  interconnection. 


Distribution  System 

The  Main  Post  distribution  system  contains  approximately  330  miles  of  pipeline  ranging  in  size  from 
6  to  30  in.,  with  smaller  individual  service  lines.  Most  of  the  lines  in  the  Main  Post  area  are  unlined  cast 
iron,  whereas  asbestos  cement  and  ductile  iron  pipe  carry  water  to  family  housing  areas.  Six  elevated  steel 
tanks  provide  a  total  system  storage  of  5  MG,  and  arc  located  at  various  points  within  the  distribution 
system.  Table  3  describes  each  tank.  Figure  3  is  a  schematic  of  the  Fort  Hood  distribution  system 
analyzed  in  this  study. 


Water  Demands 

Water  consumption  at  Fort  Hood  varies  seasonally,  with  peak  usage  occurring  during  the  summer. 
Table  4  shows  the  average  daily  water  demand  per  month  between  March  1987  and  February  1988. 
Increased  water  consumption  during  the  summer  months  is  due  in  part  to  irrigation  of  athletic  fields,  a  golf 
course,  and  lawns  in  the  family  housing  areas.  Increased  military  training  also  accounts  for  added  summer 

usage. 


Tabic  2 

Fort  Hood  Main  Pump  Station  Tank  Level  Settings* 


Pump  No. 

Level  To  Turn 

Pump  On 

Level  To  Turn 

Pump  Off 

(ft) 

(Elcv.)** 

(ft) 

(Elcv.)** 

1 

31 

1078.2 

35 

1082.2 

2 

33 

1080.2 

40 

1087.2 

3 

32 

1079.2 

39 

1086.2 

4 

22 

1069.2 

25 

1072.2 

5 

25 

1072.2 

33 

1080.2 

‘Levels  shown  represent  pumping  policy  for  summer  season. 
“Elevation  to  the  top  of  the  water  in  the  tank  in  ft  mean  sea  level  (msl). 
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Figure  2.  Fort  Hood  ground  storage  and  pumping  facility  configuration. 
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Tabic  3 


Fort  Hood  Storage  Tank  Description 


Tank  No.  Node  No.* 

Location 

Capacity 
<1000  gal) 

Max 

Elev.** 

Min 

Elev.** 

1 

25 

Pershing  Park 

500 

1086.34 

1049.00 

2 

32 

East  Post 

1000 

1086.34 

1051.59 

3 

39 

Mid  Post 

500 

1086.34 

1047.17 

4 

45 

Railhead 

1500 

1086.34 

1051.34 

5 

50 

West  Post 

500 

1086.34 

1047.17 

6 

57 

Comanche  Village 

1000 

1086.34 

1051.34 

Total  Storage:  5000 

*See  Figure  3. 
**In  ft  msl. 


Water  consumption  at  Fort  Hood  also  varies  temporally  with  increased  water  usage  during  the  midday 
and  evening  hours.  To  identify  the  diurnal  loading  pattern,  tank  water  level  variations  were  examined 
under  known  operating  conditions  (i.e.,  which  pumps  were  running).  Only  those  days  for  which  the 
optimal  control  was  applied  (30  July  1988  and  1  August  1988)  were  analyzed  for  temporal  water  variation. 
Table  5  lists  the  percentage  of  average  daily  demand  distributed  over  4-hr  periods  for  each  day  analyzed. 


Computer  Model 

The  optimal  control  methodology  combines  a  nonlinear  optimization  algorithm  with  a  simulation 
routine  in  an  attempt  to  achieve  least-cost  pumping.  The  simulation  routine,  WADISO,  solves  the 
equations  of  continuity  and  energy  associated  with  water  distribution  networks.  A  mathematical  model 
of  the  Fort  Hood  water  distribution  system  was  constructed  so  that  pump  heads  and  discharges,  tank  water 
levels,  and  pressure  at  system  nodes  could  be  determined. 


System  Schematic 

Before  a  computer  program  could  be  developed,  a  schematic  of  the  network  had  to  be  devised.  The 
resulting  schematic  illustrates  the  network  through  the  use  of  circles  and  line  segments.  Each  line  segment 
corresponds  to  an  individual  pipe  in  the  distribution  system.  Each  circle  represents  a  system  node  which 
is  either  the  location  where  two  or  more  pipes  intersect  or  where  demands  are  withdrawn  from  the  system. 
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Figure  3.  Fort  Hood  water  distribution  system  configuration. 


Table  4 


Average  Dally  Demand  per  Month 


Month/Year 

Avg.  Daily  Demand 
(gpm) 

March 

1987 

25,937 

April 

1987 

7,647 

May 

1987 

7,806 

June 

1987 

8,335 

July 

1987 

10,293 

August 

1987 

12,903 

September 

1987 

7,895 

October 

1987 

7,554 

November 

1987 

6,188 

December 

1987 

5,613 

January 

1988 

6,821 

February 

1988 

5,854 

Table  5 

Temporal  Distribution  of  Average  Daily  Demand 


Time  (hrs) 

30  July  1988 

1  Aug  1988 

Water 

Consumed  % 

Global  Demand 

Factor 

Water 

Consumed  (%) 

Global  Demand 

Factor 

0000-0400 

13.25 

1.4753 

12.71 

1.9358 

0400-0800 

16.76 

1.8653 

16.85 

2.5661 

0800-1200 

20.83 

2.3189 

20.39 

3.1062 

1200-1600 

17.39 

1.9364 

18.11 

2.7587 

1600-2000 

17.72 

1.9721 

18.93 

2.8829 

2000-2400 

14.05 

1.5642 

13.01 

1.9806 

In  computer  modeling,  it  is  usually  not  necessary  to  include  every  pipe  in  the  distribution  system. 
Instead,  it  is  possible  to  analyze  a  skeletal  system  that  includes  only  the  major  mains.  Results  of  a  skeletal 
model  arc  still  accurate  because  the  mains  not  included  in  the  skeletal  model  do  not  carry  a  great  deal  of 
flow.16  Figure  3  is  a  schematic  of  the  computer  model  developed  for  the  Fort  Hood  system. 


,#T.  M.  Walski,  Analysis  of  Water  Distribution  Systems  (Van  Nostrand  Reinhold,  1984). 
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Network  Data 


Once  a  network  schematic  of  the  system  was  developed,  a  data  base  had  to  be  established.  The  data 
base  required  information  on  each  pipe  and  node  included  in  the  schematic,  and  on  all  pumps,  tanks,  and 
flow  control  valves.  In  addition,  the  demands  associated  with  the  various  nodes  had  to  be  determined. 

Pipe  Data 

Maps  showing  the  locations  of  all  mains  in  the  system  were  obtained  from  the  Fort  Hood  Directorate 
of  Engineering  and  Housing  (DEH).  Pipes  to  be  included  in  the  model  were  selected  based  on  size  and 
importance.  Lengths  between  nodes  were  scaled  off  the  maps  and  recorded.  Roughness  values  were 
determined  by  field  tests  and  calibration  analysis.  All  pipe  data  used  in  the  analysis  are  listed  in  Appendix 
D. 

Node  Data 

Nodes  were  placed  at  locations  where  two  or  more  pipes  join  together.  Node  data  requirements 
included  elevation  and  water  use.  Elevations  for  nodes  in  the  Fort  Hood  system  were  taken  from  contour 
maps,  with  contours  spaced  at  5-ft  intervals.  Daily  water  consumption  was  estimated  from  historical 
records  provided  by  the  Fort  Hood  DEH.  Hourly  water  consumption  percentages  were  obtained  from  an 
ar.  ,ysis  of  the  operating  conditions  at  the  installation. 

In  actual  operation,  water  is  withdrawn  from  a  water  main  at  various  service  connections  located 
along  the  entire  length  of  the  pipe.  However,  when  modeling  a  distribution  system,  these  demands  are 
aggregated  at  individual  nodes  in  the  immediate  vicinity  of  water  withdrawal.  For  example,  the  combined 
water  demand  for  an  entire  family  housing  area  may  be  placed  at  a  single  node. 

In  modeling  the  Fort  Hood  distribution  system,  the  average  daily  demand  was  distributed  among 
several  nodes.  The  selection  of  demand  nodes  and  the  amount  of  water  distributed  to  each  node  were 
based  on  an  analysis  of  major  users  at  the  Main  Post.  These  users  included  family  housing,  barracks, 
washracks,  the  east  post  airfield,  and  Damall  Hospital.  A  complete  set  of  node  data  is  provided  in 
Appendix  D. 

Tank  Data 

Six  storage  tanks  were  included  in  the  computer  model.  All  tanks  operate  in  the  same  pressure  zone 
and  have  the  same  overflow  elevation  (Table  3). 

PRV  Data 

At  present,  the  Fort  Hood  system  contains  two  pressure-reducing  valves  (PRVs)  isolating  the  Pershing 
Park  and  Venable  Village  areas  (nodes  19  through  23  in  Figure  3)  from  the  rest  of  the  system.  The  PRVs 
arc  necessary  to  reduce  the  pressures  in  these  areas  since  they  lie  at  low  elevations  with  respect  to  the  rest 
of  the  system  and  are  close  to  the  main  pump  station.  The  hydraulic  grade  of  both  PRVs  was  set  at  1 101 
ft. 


Pump  Field  Tests 

As  part  of  the  optimal  control  study,  the  pumps  at  the  main  pump  station  were  field-tested  to 
determine  their  actual  operating  characteristics.  The  testing  was  conducted  on  17  and  18  March  1988. 
Each  pump  was  tested  individually  and,  in  some  cases,  parallel  combinations  of  pumps  were  tested  to 
obtain  data  over  a  wide  variety  of  operating  conditions. 
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The  main  pump  station  contains  live  vertical  turbine  pumps,  liacli  pump  is  rated  at  202  f  t  of  head 
(Table  1),  but  the  discharges  vary.  All  pumps  pull  water  from  a  common  suction  line  and  discharge  water 
into  a  single  discharge  header. 

Testing  Procedure 

For  each  test,  the  pump  was  started  against  a  closed  valve  in  the  discharge  line  to  obtain  a  shutoff 
head.  After  the  pump  reached  full  speed,  a  pressure  reading  was  taken  from  a  gauge  located  in  the 
discharge  line  between  the  pump  and  the  valve.  A  pressure  gauge  installed  at  the  base  of  the  old  500-hp 
pump  provided  a  reading  for  suction  head.  The  difference  in  pressure  between  the  two  gauge  displays 
was  the  head  delivered  by  the  pump.  The  power  consumption  of  the  pump  was  also  recorded.  Once 
shutoff  head  and  corresponding  power  were  obtained,  the  valve  was  opened  slightly  and  pressure, 
discharge,  and  power  consumption  values  were  recorded.  This  process  was  repeated,  with  the  valve 
opening  increased  each  time  until  the  valve  was  completely  open. 

Pump  discharge  was  measured  using  a  pitot  tube  inserted  in  the  pump  discharge  header.  The  pitot 
tube  provided  a  pressure  differential  that  was  converted  into  pipeline  velocity  and  ultimately  into  pipe  flow 
rate.  Power  readings  in  the  form  of  current  draw  and  voltage  drop  across  the  pump  were  taken  and 
converted  into  kilowatt  consumption. 

Pump  Data 

For  this  analysis,  only  the  pumps  located  at  the  main  pump  station  were  included.  Although  none  of 
the  booster  pumps  located  at  Fort  Hood  were  considered  in  the  study,  the  proposed  methodology  can 
easily  handle  such  situations.  To  include  booster  stations,  the  pump’s  characteristics  must  be  known  and 
they  must  use  electrical  power. 

A  careful  review  of  the  data  indicated  that  nearly  all  of  the  pumps  have  characteristics  close  to  the 
manufacturers’  original  specifications.  As  a  result,  the  pump  data  provided  by  the  manufacturers  were 
used  in  the  analysis.  The  old  500-hp  pump  had  pump  heads  and  discharges  different  from  those  shown 
by  the  original  curve.  Consequently,  this  pump  was  omitted  from  the  analysis. 


System  Calibration 

When  using  a  computer  simulation  routine  to  analyze  a  water  distribution  system,  it  is  extremely 
important  that  the  mathematical  model  of  the  system  be  an  accurate  representation  of  actual  field 
conditions.  If  this  is  not  the  case,  the  results  provided  by  the  simulation  routine  will  be  of  limited  value. 
Therefore,  the  model  must  be  calibrated.  Calibration  is  done  by  adjusting  both  water  usage  and  pipe 
roughness  until  heads  and  flows  computed  by  the  simulation  routine  match  those  observed  in  field 
measurements. 

One  common  method  of  obtaining  calibration  data  is  to  conduct  fire  hydrant  flow  tests  at  various 
locations  throughout  the  system.  These  tests  provide  hydraulic  information  (pressures  and  flows)  during 
both  low  flow  and  high  flow  periods.  This  information  can  be  used  with  the  network’s  computer  model 
to  calibrate  the  system. 

Calibration  data  for  the  Fort  Hood  Main  Post  distribution  system  were  gathered  from  fire  hydrant  flow 
tests  conducted  on  15  and  16  March  1988.  The  location  of  each  fire  hydrant  flow  test  is  shown  in  Figure 
4;  Tabic  6  lists  results  for  each  test.  For  each  test,  at  least  two  and,  in  most  cases,  three  hydrants  were 
flowed  to  induce  as  much  stress  into  the  system  as  possible. 
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Figure  4.  Locations  of  fire  hydrant  flow  tests  conducted  at  Fort  Hood. 


Table  6 


Hydrant  Flow  Test  Results  for  Fort  Hood  Main  Post 


Combined 

Number  of 

Hydrant 

Residual 

Kiev. 

Hydrants 

Pressure 

Discharge 

Test  No. 

Location 

Node 

(ft) 

Flowed 

(psi) 

(«pm) 

1 

Pershing  Park 

70 

878 

0 

89.5 

0 

1 

61.5 

1288 

2 

40 

1922 

3 

28 

2386 

2 

Comanche  Village 

80 

888 

0 

82 

0 

1 

77 

1336 

2 

70 

2357 

3 

62 

3320 

3 

Duncan  School 

81 

923 

0 

68 

0 

1 

58 

1033 

2 

38 

1708 

4 

Anderson  Golf 

27 

905 

0 

83 

0 

Course 

1 

54 

903 

2 

30 

1254 

5 

Meadows  School 

72 

906 

0 

83 

0 

1 

77 

1277 

2 

71 

2509 

3 

70 

3717 

6 

East  Fort  Hood 

35 

891 

0 

84 

0 

Airfield 

1 

75 

1351 

2 

72 

2120 

3 

68 

3283 

7 

McNair  Village 

75 

912 

0 

84 

0 

1 

74 

1220 

2 

71 

2451 

3 

69 

3384 

8 

Central  Avc. 

74  ■ 

924 

0 

69 

0 

between  37th 

1 

63 

1185 

and  42nd  St. 

2 

57 

2112 

2 

57 

2112 

9 

Battalion  Avc 

76 

937 

0 

62 

0 

between  52nd  and 

1 

62 

1060 

62nd  St. 

2 

61 

2184 

3 

58 

3318 

10 

Warehouse  Rd. 

79 

965 

0 

49 

0 

and  72nd  St. 

1 

46 

1073 

2 

41 

1936 

3 

37 

2701 

11 

Support  Avc.  and 

78 

930 

0 

66 

0 

72nd  St. 

1 

65 

1099 

2 

64 

2350 

3 

63 

3490 

12 

Battalion  Avc. 

55 

898 

0 

80 

0 

and  78lh  St. 

1 

75 

1320 

2 

67 

2542 

3 

60 

3495 
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The  WADISO  model  of  the  water  distribution  system  was  calibrated  using  a  nonlinear  optimization 
technique.17  The  underlying  theory  of  this  technique  is  very  similar  to  the  optimal  control  methodology 
used  to  develop  the  system.  The  optimization  code,  GRG2  was  coupled  with  a  hydraulic  simulation 
routine,  WADISO,  to  minimize  the  difference  between  observed  and  predicted  system  heads  and  observed 
and  predicted  tank  water  levels.  The  decision  variables  in  the  model  are  pipe  roughness,  nodal  water 
demand,  and  global  demand  factor. 

Pipe  roughness  and  nodal  water  demand  are  systematically  adjusted  within  defined  limits  until  the 
difference  between  observed  and  predicted  heads  is  minimized.  The  global  demand  factor  is  a  constant 
by  which  each  nodal  demand  is  multiplied  to  account  for  the  hourly  variation  in  water  demand.  This 
variable  is  also  adjusted  between  defined  limits  until  the  difference  between  observed  and  predicted  tank 
water  levels  is  minimized. 

The  calibration  technique  allows  the  adjusted  parameters  to  be  constrained  between  upper  and  lower 
values  which  are  selected  based  on  a  knowledge  of  the  system.  Table  7  compares  the  observed  and 
predicted  system  heads  after  calibration  for  selected  fire  flow  tests.  Calibration  results  for  the  Fort  Hood 
system  are  presented  in  Appendix  D  as  Hazen-Williams  C-Factor  and  Node  Demand.  Calibrated  global 
demand  factors  are  shown  in  Table  5. 

Note  that  only  one  tank  in  the  Fort  Hood  distribution  system,  tank  5,  is  monitored  for  variation  in 
water  level.  Thus,  information  on  the  water  levels  for  the  remaining  tanks  was  not  available.  For  the 
global  demand  factor  calibration,  tanks  1  through  4  were  assumed  to  be  full  at  the  beginning  of  the 
analysis,  the  level  in  tank  5  was  known,  and  the  level  in  tank  6  was  assumed  to  be  5  ft  less  than  the  level 
in  tank  5.  These  assumed  starting  elevations  were  chosen  based  on  the  results  of  the  stcadystatc 
simulation  of  the  Fort  Hood  system  and  on  conversations  with  installation  personnel. 

The  actual  lank  levels  over  the  course  of  the  day  for  all  tanks  other  than  tank  5  were  obtained  from 
the  global  demand  factor  calibration.  The  starting  elevations  for  all  tanks  other  than  tank  5  were  assumed. 
The  global  demand  factor  for  a  given  time  period  was  found  by  minimizing  the  difference  between  the 
known  and  predicted  levels  for  tank  5.  The  ending  elevation  of  each  tank  was  recorded  and  used  as  the 
beginning  level  for  the  next  time  step.  This  procedure  was  repeated  for  each  time  step  analyzed. 


r  W.  W.  Sharp  and  D.  W.  Chase,  "Verifying  Water  Distribution  Computer  Model  Calibration  Results,"  /WWV/t  Computers  and 
Automation  Specialty  Conference ,  Denver,  CO  (1989). 
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Tabic  7 


Calibration  Results  for  Fort  Hood  Main  Post 


Location 

Number  of 

Hydrants  Flowed 

Observed 

Head 

(ft) 

Predicted 

Head 

(ft) 

1 

Pershing  Park 

0 

1085 

1088 

2 

970 

970 

2 

Comanche  Village 

0 

1077 

1079 

1 

1066 

1057 

2 

1050 

1047 

3 

1031 

1030 

3 

Duncan  School 

0 

1080 

1079 

1 

1057 

1048 

2 

1011 

1007 

4 

Meadows  School 

0 

1098 

1101 

1 

1084 

1087 

3 

1068 

1068 

5 

East  Fort  Hood 

0 

1085 

1097 

Airfield 

1 

1064 

1077 

2 

1057 

1065 

1048 

1041 

6 

McNair  Village 

1 

1083 

1090 

2 

1076 

1087 

7 

Central  Ave.  between 

0 

1083 

1092 

37th  and  42nd  St. 

1 

1070 

1075 

2 

1056 

1055 

8 

Battalion  Ave.  between 

0 

1080 

1089 

52nd  and  62nd  St. 

1 

1080 

1084 

2 

1079 

1081 

3 

1071 

1080 

9 

Warehouse  Rd.  and  72nd  St. 

0 

1078 

1080 

1 

1071 

1074 

2 

1060 

1064 

3 

1050 

1051 

10 

Support  Ave.  and  72nd  St. 

0 

1082 

1079 

1 

1080 

1077 

2 

1078 

1077 

3 

1076 

1076 

1  1 

Battalion  Ave.  and 

0 

1083 

1078 

78th  St. 

1 

1071 

1067 
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4  APPLICATION  OF  THE  OPTIMAL  CONTROL  METHODOLOGY 


The  optimal  control  methodology  was  applied  to  the  Fort  Hood  water  distribution  system  to  determine 
the  effectiveness  of  the  model  and  demonstrate  its  capabilities.  The  methodology  was  applied  to  2  days 
for  which  data  on  actual  operating  conditions  were  available.  The  days  were  30  July  1988  (weekend)  and 
1  August  1988  (weekday). 

A  review  of  actual  lank  levels  and  pump  operating  times  indicated  that  Fort  Hood  system  demands 
were  cyclic  with  a  period  of  1  day.  The  demands  were  approximately  the  same  during  the  week  from 
Monday  to  Friday.  System  demands  were  less  on  the  weekends;  however,  the  weekend  demands  were 
also  cyclic  with  a  daily  period.  For  this  reason,  2  days,  one  during  the  week  and  the  other  on  the 
weekend,  were  selected.  A  review  of  climatic  data  indicated  that  precipitation  and  temperature  were 
normal  for  several  days  up  to  and  including  the  days  analyzed. 


Electrical  Rate  Structure 

Fort  Hood  purchases  electricity  from  Texas  Utilities  Electric  Company  under  a  general  service  rate 
schedule  at  an  average  cost  of  $0.045/kWh.  The  contract  oetween  Fort  Hood  and  the  Texas  Utilities 
Electric  Company  contains  a  time-of-day  option;  however,  this  option  applies  only  to  electrical  demand 
charges  and  not  to  energy  unit  charges.  For  this  study,  the  average  cost  of  electricity  is  applied  evenly 
over  the  course  of  the  day. 


Actual  Operating  Results 

Pump  operation  at  Fort  Hood  is  controlled  by  the  water  level  in  tank  5  or  the  West  Post  Tank. 
Therefore,  this  is  the  only  tank  in  which  the  water  levels  are  monitored.  Analysis  of  the  calibration  results 
indicated  that  tanks  1  and  4  stay  completely  full  during  normal  operating  conditions.  This  was  verified 
from  conversations  with  Fort  Hood  personnel.  During  periods  of  high  demand  or  when  the  100-hp  pump 
is  operating  independently,  tanks  1  and  4  will  drain.  Because  of  its  location,  water  must  travel  through 
many  feet  of  smatl  diameter  pipe  between  the  pump  station  and  tank  6.  The  resulting  head  loss  causes 
tank  6  to  drain  under  normal  operating  conditions. 

Actual  Cost  of  Operations 

The  Fort  Hood  pump  station  does  not  have  electric  meters  installed  on  the  pumps.  As  a  result,  the 
actual  electrical  consumption  and  hence  the  operating  cost  cannot  be  computed  directly.  The  actual 
operating  cost  used  in  this  study  was  computed  based  on  the  known  pump  operating  times,  the  known  tank 
trajectory  for  tank  5,  and  the  calibrated  tank  trajectories  for  the  remaining  tanks.  As  a  result,  the  cost  is 
not  a  true  cost  of  operation,  but  is  a  computed  cost  based  on  available  information. 

30  July  I9SS 

The  water  level  in  tank  5  at  midnight  was  36.12  ft  and  dropped  to  a  low  level  of  31  ft  in  the  early 
morning  and  late  afternoon  hours.  The  water  level  at  the  end  of  the  day  was  32  ft.  During  the  day,  the 
100-.  200- .  and  300-hp  pumps  were  operated.  Table  8  lists  the  pump  operating  times  and  corresponding 
tank  water  levels  for  tank  5.  Note  that  pump  3  cycles  on  and  off  approximately  every  hour  during  the 
early  morning  hours.  The  actual  cost  of  operation  for  Saturday,  30  July  1988,  was  $380.36. 
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Table  8 


Actual  I’ump  Operating  Times  for  30  July  1988 


Time 

(hr) 

Beginning 

Elevation 

(ft) 

Ending 

Elevation 

(ft) 

Pumps 

Operating 

0000-0030 

36.12 

32.00 

2 

0030-0145 

32.00 

39.00 

2,3 

0145-0230 

39.00 

32.00 

2 

0230-0345 

32.00 

39.00 

2,3 

0345  0445 

39.00 

32.00 

2 

0445-0900 

32.00 

31.00 

2,3 

0900-1500 

31.00 

35.00 

1.2.3 

1500-1830 

35.00 

31.00 

2,3 

1830  21 20 

21.00 

35.00 

1,2,3 

2120  2300 

35.00 

39.00 

2.3 

2300-2400 

39.00 

32.00 

2 

/  August  19XX 

The  beginning  and  ending  water  levels  for  tank  5  were  35.18  and  35.00  ft,  respectively.  The  low 
water  level  for  the  day  was  25  ft.  The  100-,  200-,  300-,  and  new  500-hp  pumps  were  operated  at  various 
times  during  the  day.  Table  9  lists  the  pump  operating  times  and  associated  tank  levels  for  tank  5  on 
Monday,  1  August  1988.  As  is  the  case  on  30  July  1988,  pumps  (3  and  5  in  this  case)  cycle  on  and  off 
several  times  during  the  day.  The  actual  cost  of  operation  for  this  day  was  $560.75. 


Optimal  Operating  Results 

OCOPS  was  applied  to  the  2  test  days  to  determine  if  the  actual  operating  procedure  at  Fort  Hood 
could  be  modified  to  achieve  a  cost  of  pumping  less  than  the  actual  cost  while  still  maintaining  sufficient 
system  pressure  and  adequate  system  storage.  The  time  interval  selected  for  analysis  was  4  hr. 
Constraints  were  placed  on  the  system’s  pressures  and  storage.  These  constraints  were  evaluated  at  the 
end  of  each  time  interval. 

Implicit  Bound  Constraints 

implicit  bound  constraints  (Appendix  A)  were  placed  on  system  pressure  and  storage.  The  pressure 
constraint  was  applied  to  nodes  24,  42,  and  47.  These  nodes  were  selected  based  on  their  elevation  and 
location  within  the  distribution  system.  Analysis  of  the  calibration  results  indicated  that  node  24  had  the 
highest  pressure  (large  diameter  main  close  to  the  pump  station)  within  the  system.  Nodes  42  and  47  had 
the  lowest  pressures  (high  elevation)  in  the  system. 

Nodes  at  the  suction  and  discharge  sides  of  the  pumps  were  not  included  in  the  pressure  constraint 
set  since  they  were  expected  to  have  very  low  and  high  pressures,  respectively.  The  maximum  allowable 
pressure  was  set  at  100  psi  while  the  minimum  acceptable  pressure  was  40  psi. 
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Table  9 


Actual  Pump  Operating  Times  for  1  August  1988 


Time 

(hr) 

Beginning 

Elevation 

(ft) 

Ending 

Elevation 

(ft) 

Pumps 

Operating 

0000-0309 

35.18 

39.00 

2.3 

0309-0340 

39.00 

32.00 

2 

0340-0600 

32.00 

31.00 

2.3 

0600-0700 

31.00 

25.00 

1.2,3 

0700-0830 

25.00 

33.00 

1. 2,3,5 

0830-0930 

33.00 

25.00 

1.2,3 

0930-1115 

25.00 

33.00 

1. 2.3,5 

1115-1200 

33.00 

25.00 

1.2.3 

1200-1345 

25.00 

32.00 

1. 2.3,5 

1345-1510 

33.00 

25.00 

1.2,3 

1510-1615 

25.00 

33.00 

1,2, 3, 5 

1615-1725 

33.00 

25.00 

1,2.3 

1725-1900 

25.00 

33.00 

1.2,3, 5 

1900-2000 

33.00 

25.00 

1,2.3 

2000-2100 

25.00 

33.00 

1. 2.3.5 

2100-2400 

33.00 

35.00 

1.2,3 

Constraints  were  also  placed  on  all  tanks  in  order  to  keep  a  reserve  of  water  in  the  system  for  fire¬ 
fighting  needs  and  to  have  the  tank  levels  at  a  desired  ending  elevation.  The  upper  bound  placed  on  the 
water  level  was  the  tank’s  top  elevation.  This  point  was  selected  so  that  the  tank  would  not  fill  past  its 
overflow  elevation.  The  upper  bound  of  the  tank  constraint  equal  to  the  top  elevation  of  the  tank  was  in 
force  for  each  time  interval. 

The  lower  bound  tank  constraint  was  dependent  on  two  items:  (1)  the  amount  of  storage  necessary 
for  fire-fighting  purposes  and  (2)  the  actual  ending  elevation  of  the  tank.  For  all  time  intervals  except  the 
last  one,  the  lower  tank  constraint  depended  on  the  storage  necessary  in  the  system.  The  available  storage 
in  the  system  during  these  time  intervals  was  constrained  to  be  no  less  than  1.25  MG.  This  amount  is 
enough  to  combat  a  5000-gpm  fire  for  4  hr.  In  most  cases,  the  available  storage  was  greater  than  1.25 
MG.  Table  10  lists  the  upper  and  lower  constraints  placed  on  each  tank. 

For  the  last  time  interval,  the  lower  bound  tank  constraint  was  set  at  the  actual  ending  elevation  of 
the  tanks  for  the  days  analyzed.  The  optimization  algorithm  attempts  to  force  tank  levels  to  their  lower 
bound  during  the  final  time  interval  since  a  lower  tank  level  results  in  a  lower  pumping  head,  which 
results  in  lower  operating  costs.  Thus,  it  is  expected  that  the  final  tank  levels  from  the  optimal  procedure 
would  be  equal  or  slightly  above  the  final  tank  levels  resulting  from  the  actual  procedure.  This  convention 
enables  the  actual  and  optimal  operating  policies  to  be  compared  equally. 

A  4-hr  time  interval  was  selected  to  reduce  the  computer  run  times.  The  disadvantage  to  using  larger 
time  increments  is  that  the  method’s  flexibility  is  compromised.  Pumps  are  actually  operated  over  a 
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continuous  timeframe  as  opposed  to  a  discrete  timeframe.  The  sinallcr  the  time  interval,  the  closer  the 
results  approach  continuous  operation.  Consequently,  the  smaller  the  time  step,  the  greater  the  possibility 
for  savings.  Unfortunately,  smaller  time  steps  increase  the  time  necessary  to  reach  an  optimal  solution 
and  could  result  in  pumps  cycling  on  and  off  frequently.  Pump  cycling  is  undesirable  since  it  wastes 
energy  and  causes  wear  and  tear  on  a  pump. 

30  July  1988 

The  optimization  algorithm  used  in  the  model  requires  an  initial  guess  to  start  the  optimization 
process.  The  guess  selected  for  this  day  was  chosen  to  approximate  the  actual  pump  operating  times. 
This  choice  resulted  in  pump  2  operating  for  the  entire  24-hr  period  and  pump  3  running  for  all  but  4.5 
hr.  Pump  1  ran  for  about  8  hr  while  pump  5  did  not  operate  at  all.  The  cost  of  operation  under  the  initial 
guess  was  $366.92.  Figure  5  compares  the  actual  operating  times  and  the  operating  times  corresponding 
to  the  initial  guess.  Note  that,  in  the  optimal  control  methodology,  pumps  are  not  allowed  to  begin 
operating  in  the  middle  of  a  time  interval.  Rather  they  must  start  operating  at  the  beginning  of  the  time 
interval  and  run  continuously  for  X, §i  hours. 

Although  the  initial  guess  resulted  in  an  operating  cost  less  than  the  actual  cost,  the  ending  tank 
levels  for  each  tank  were  less  than  the  ending  levels  resulting  from  the  calibration  analysis.  No  pressure 
constraints  were  violated,  however.  The  ending  tank  level  constraint  violation  can  be  attributed  to  the 
initial  guess  only  approximating  actual  pump  operating  times  and  to  the  global  demand  factor  being  a  4-hr 
average  computed  from  calibrated  global  demand  factors  in  effect  throughout  the  day. 

The  optimization  algorithm  entered  into  a  Phase  1  optimization  to  satisfy  the  violated  tank 
constraints.  The  constraints  were  satisfied  by  increasing  the  pump  run  times  for  all  pumps  not  already 
at  a  lower  or  upper  bound.  The  cost  of  pumping  after  all  constraints  were  satisfied  was  $392.02.  Once 
all  constraints  were  satisfied,  the  algorithm  entered  into  a  Phase  2  optimization  to  complete  the 
optimization. 

The  cost  of  operation  provided  by  the  OCOPS  upon  completion  of  the  Phase  2  optimization  was 
$390.00.  This  cost  is  approximately  $10  greater  than  the  actual  cost  of  operation.  No  constraints  were 
violated  under  the  final  solution.  Most  of  the  ending  elevatioas  for  the  tanks  were  close  to  their  lower 


Tabic  10 

Tank  Water  Level  Constraints 


Minimum  Acceptable 

Maximum  Allowable 

Water  Level 

Water  Level 

Tank  No. 

(rt) 

(ft) 

1 

9.00 

37.34 

2 

9.00 

34.75 

3 

9.00 

39.17 

4 

9.00 

35.00 

5 

9.00 

39.17 

6 

9.00 

35.00 
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ACTUAL  PUMP  OPERATING  TIMES 
FOR  30  JULY  1988 
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INITIAL  OPTIMIZATION  GUESS 
FOR  30  JULY  1988 
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Figure  5.  Actual  and  initial  guess  operating  times  for  30  July  1988. 
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bound.  This  result  was  expected  because  the  nature  of  the  problem  causes  tank  levels  to  seek  their  lower 
bounds  since  pumping  heads  are  a  function  of  tank  water  level.  The  lower  the  tank  level,  the  lower  the 
pump  head,  and  hence  the  lower  the  cost  of  operation.  Figure  6  depicts  the  pump  operating  times  for  the 
final  solution,  and  Figures  7  through  12  show  the  resulting  tank  trajectories  for  actual  and  optimal 
conditions. 

1  August  1988 

The  initial  guess  provided  for  this  day  again  was  an  approximation  of  the  actual  operating  times 
for  pumps.  This  guess  resulted  in  pumps  2  and  3  operating  for  the  entire  day,  pump  1  operating  for  18 
hr,  and  pump  5  running  for  about  8  hr.  The  cost  of  operation  associated  with  the  initial  guess  was 
$555.82.  Figure  13  compares  the  actual  operating  times  and  the  time  associated  with  the  initial  guess  for 
1  August  1988. 

As  was  the  case  with  30  July  1988,  the  initial  guess  resulted  in  the  final  water  level  in  several 
tanks  being  below  the  desired  ending  elevation.  As  a  result,  the  optimization  algorithm  entered  into  a 
Phase  1  optimization  to  satisfy  the  tank  level  constraints.  This  was  accomplished  by  increasing  the  value 
of  decision  variables  not  already  at  a  bound.  No  other  pressure  or  tank  level  constraints  were  violated  due 
to  the  initial  guess. 

The  cost  of  pump  operation  provided  by  the  optimal  solution  was  $562.88  which  was  about  $2 
greater  than  the  actual  cost  of  pumping.  Like  the  previous  case,  all  constraints  were  satisfied  and  the 
ending  tank  water  level  for  all  tanks  was  either  at  or  slightly  above  its  lower  bound.  Figure  14  shows  the 
operating  times  for  each  pump  corresponding  to  the  optimal  solution.  Figures  15  through  20  show  the 
tank  levels  associated  with  actual  and  optimal  operation. 


Discussion 

The  benefit  of  applying  mathematical  optimization  to  pump  operating  strategies  for  water  distribu¬ 
tion  systems  has  been  demonstrated  previously.18  For  example,  an  optimization  approach  using  dynamic 
programming  (DP)  was  applied  to  the  Washington,  DC  system  and  resulted  in  projected  annual  savings 
of  nearly  $100,000.19  DP  works  quite  well  for  systems  such  as  that  in  Washington  which  have  only  one 
or  two  storage  tanks  in  a  service  area.  However,  most  systems  have  more  than  one  or  two  tanks  in  a 
service  area  and  as  a  result,  the  use  of  DP  becomes  infeasible  due  to  dimensionality.  In  this  study,  the 
nonlinear  programming  algorithm  OCOPS  was  developed  for  use  with  highly  dimensional  systems  such 
as  the  one  at  Fort  Hood,  TX. 

All  optimization  models  require  accurate  simulation  of  system  response.  In  the  case  of  pump 
optimization,  the  water  distribution  system  response  is  described  in  the  form  of  the  pressure  and  flow 
variation  within  the  network,  the  changes  in  tank  level,  and  values  for  pump  head  and  pump  discharge. 
The  modeled  response  must  be  very  close  to  the  actual  system  response  for  the  optimization  algorithm  to 
operate  correctly  and  the  simulation  routine  to  provide  meaningful  results.  In  other  words,  the  simulation 
model  must  be  well  calibrated  for  the  system  being  optimized. 

As  discussed  earlier,  the  lack  of  information  on  tank  level  variation  prevented  complete  calibration 
of  the  Fort  Hood  system.  These  data  were  available  for  one  tank  only  at  Fort  Hood.  Tank  trajectories 
for  the  remaining  tanks  were  not  known.  Therefore,  the  starting  elevation  for  five  of  the  six  tanks  had 
to  be  assumed  and  the  system  calibrated  based  on  information  for  the  one  tank. 


llfW.  H.  Clingcnpeel;  A.J.  Tarquin  and  J.  Dowdy. 
'*L.  Ormsbee. 
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Figure  11.  Tank  5  water  levels  for  30  July  1988. 
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Figure  12.  Tank  6  water  levels  for  30  July  1988. 
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ACTUAL  PUMP  OPERATING  TIMES 
FOR  1  AUG  1988 


Figure  13.  Actual  and  initial  guess  operating  times  for  1  August  1988. 
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Figure  14.  Optimal  pump  operating  times  for  1  August  1988, 
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Figure  15.  Tank  1  water  levels  for  1  August  1988. 
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Figure  16.  Tank  2  water  levels  for  1  August  1988. 
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Figure  17.  Tank  3  water  levels  for  1  August  1988. 


Although  the  OCOPS  results  did  not  show  an  improvement  over  the  estimated  actual  cost  of 
operation,  some  comments  can  be  made  on  the  model  as  applied  to  Fort  Hood.  Several  initial  guesses  for 
pump  operating  time  were  used  during  the  study.  These  starting  points  ranged  from  ones  with  all  pumps 
running  for  the  entire  day  to  those  with  only  two  pumps  running  for  a  portion  of  the  day.  In  all  cases, 
OCOPS  converged  to  a  solution  within  a  few  dollars-well  within  the  accuracy  of  the  simulalion-of  the 
estimated  actual  cost.  This  result  is  encouraging  since  the  model  was  able  to  reach  the  same  solution, 
regardless  of  the  initial  guess  supplied. 

Several  reasons  could  explain  why  improvements  to  the  estimated  actual  cost  were  not  shown  for  the 
Fort  Hood  system.  The  use  of  a  questionably  calibrated  mathematical  model  of  the  Fort  Hood  system 
almost  certainly  influenced  the  results.  Another  reason  could  be  that  the  Fort  Hood  system  may  be  so 
highly  constrained  and  have  such  a  narrow  feasible  path  that  the  solution  has  been  reached  heuristically 
by  trial  and  error.  Although  the  system  has  5  MG  of  elevated  storage,  it  is  distributed  throughout  the 
system.  There  may  not  be  enough  available  storage  in  any  one  tank  to  provide  the  system  with  the 
flexibility  needed  for  improved  pump  operation.  Rather,  the  storage  may  be  needed  to  supply  system 
demands  immediately  instead  of  storing  water  for  future  use.  Finally,  the  Fort  Hood  system  may  be  veiy 
sensitive  to  the  time  pumps  are  operating  If  so,  it  may  have  been  more  realistic  to  use  a  smaller  time 
interval  than  the  4-hr  time  step  used  in  GiOer  to  improve  actual  operation.  Reducing  the  time  interval 
would  have  reduced  the  constraint  placed  on  time  and  provided  flexibility  by  more  closely  approximating 
continuous  operation. 


Summary 

The  results  from  OCOPS  were  initially  surprising  since  the  cost  provided  by  the  model  was 
expected  to  be  less  than  the  estimated  actual  cost  of  operation.  However,  due  to  the  questionable 
calibration,  it  is  unfair  to  speculate  on  the  model’s  capabilities  without  first  testing  it  on  a  well  calibrated 
system  for  which  real  actual  costs  are  known.  It  is  possible  that  the  Fort  Hood  system  has  been  optimized 
heuristically.  If  so,  then  the  results  of  OCOPS  are  indeed  encouraging  since  they  are  very  close  to  the 
"optimal'1  results.  Before  extolling  the  virtues  of  the  model,  however,  OCOPS  should  be  tested  using  a 
system  with  a  truly  calibrated  simulation  model.  Th„  results  for  the  Fort  Hood  system  must  be  treated 
with  caution  since  the  calibration  of  the  water  distribution  system  was  questionable.  To  use  OCOPS  or 
any  pump  optimization  model  with  confidence,  the  model  must  be  able  to  simulate  the  response  of  the 
water  distribution  system  accurately.  If  the  modeled  response  of  the  system  is  not  accurate,  then  the 
results  of  the  optimization  model  may  be  misleading. 

Although  the  results  of  this  study  indicate  that  costs  provided  by  OCOPS  are  no  less  than  the 
estimated  actual  cost  of  operation  at  Fort  Hood,  the  technique  appears  to  work.  It  is  possible  that  room 
for  improvement  exists  at  Fort  Hood,  but  without  a  well  calibrated  model,  the  magnitude  of  improvement 
is  impossible  to  quantify. 


Advantages  of  Pump  Optimization  Models 

In  general,  water  distribution  systems  are  highly  dimensional  with  multiple  storage  tanks  and 
pressure  zones.  Most  water  systems  have  adequate  storage  to  meet  peak  demands  and  provide  reserve 
water  for  emergencies  or  fire  fighting.  However,  in  some  cases,  these  same  systems  may  be  highly 
constrained.  The  storage  available  in  these  systems  may  be  insufficient  to  store  elevated  water  for  use  at 
later  times. 

A  good  optimization  model  capable  of  handling  multiple  storage  tanks,  coupled  with  a  well 
calibrated  model  of  the  water  distribution  system,  can  identify  the  potential  for  improved  operation.  Such 
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improved  operation  cannot  only  surface  in  the  form  of  reduced  electrical  costs,  but  also  in  more  efficient 
system  operation  since  these  models  generally  have  reliability  constraints  embedded  within  them. 

It  should  be  emphasized  that,  for  some  systems,  a  small  percentage  reduction  in  energy  costs  can 
translate  into  a  large  dollar  savings.  A  pump  optimization  model  can  identify  pump  combinations  that 
result  in  minimal  energy  costs.  It  is  possible  that  some  systems  may  be  operating  optimally  without  the 
use  of  a  model.  However,  in  these  cases,  an  optimization  model  can  still  provide  important  information 
on  how  the  system  reacts  to  changes  in  system  operation,  storage,  demand,  electrical  rate  structure,  and 
piping  configuration. 
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5  CONCLUSIONS  AND  RECOMMENDATIONS 


An  optimization  methodology  to  minimize  pumping  costs  has  been  developed  and  applied  to  the  Fort 
Hood,  TX  water  distribution  system.  The  methodology  combines  a  nonlinear  optimization  algorithm, 
GRG2,  with  a  water  distribution  simulation  routine,  WADISO.  The  two  programs  are  linked  together  by 
a  GCGiVlP  routine  which  computes  the  values  of  die  objective  function  and  the  problem  constraints.  The 
combined  subroutines,  along  with  a  data  entry  routine,  are  contained  in  a  computer  program  titled  Optimal 
Control  of  Pump  Stations  (OCOPS). 

OCOPS  adjusts  the  fraction  of  time  a  pump  operates  during  a  particular  interval  to  minimize  electrical 
costs  while  maintaining  system  reliability.  System  reliability  is  described  in  the  form  of  problem 
constraints  on  system  storage  and  pressure.  Additional  constraints  on  the  problem  include  explicit  bound 
constraints  on  decision  variables  and  implicit  system  constraints. 

Implicit  system  constraints  are  the  equations  describing  the  head  and  flow  distribution  in  a  hydraulic 
network  and  are  solved  by  the  WADISO  program.  The  remaining  constraints  are  solved  by  GRG2.  In 
OCOPS,  the  implicit  bound  constraints  on  system  parameters  include  an  upper  and  lower  bound  on  nodal 
pressures  and  a  maximum  and  minimum  allowable  tank  level.  Appendix  A  provides  a  detailed  description 
of  the  problem  formulation. 

The  optimization  technique  was  applied  to  the  Fort  Hood  water  distribution  system  to  test  its 
capabilities.  The  model  was  applied  to  2  days  for  which  actual  operating  conditions  were  known:  30  July 
1988  and  1  August  1988.  The  actual  cost  of  pumping  for  30  July  1988  was  $380.36,  whereas  the  cost 
of  pumping  provided  by  the  model  was  $390.00.  For  1  August  1988,  the  actual  power  cost  was  $560.75 
and  the  cost  provided  by  the  model  was  $562.88.  All  problem  constraints  were  satisfied  under  the  final 
solution. 

OCOPS  was  unable  to  improve  the  current  pumping  schedule  for  the  Fort  Hood  system.  Most  likely 
this  is  due  to  a  poorly  calibrated  model  of  the  Fort  Hood  system.  It  is  possible,  however,  that  the  Fort 
Hood  system  is  so  highly  constrained  with  respect  to  system  storage  that  it  may  already  be  operating 
optimally. 

Even  though  OCOPS  was  unable  to  provide  a  pump  schedule  for  Fort  Hood  which  resulted  in 
operating  costs  less  than  the  estimated  actual  costs,  the  model  should  not  be  abandoned.  Rather,  the 
program  should  be  tested  on  a  system  for  which  a  well  calibrated  mathematical  model  is  available  to  test 
the  model’s  true  capabilities.  The  model  docs  appear  to  work  since  it  converged  to  a  solution  very  close 
to  the  estimated  actual  cost  for  a  variety  of  initial  pump  schedules  supplied. 

A  useful  feature  of  OCOPS  is  its  ability  to  handle  systems  with  multiple  storage  tanks.  Previous 
methods  have  generally  relied  on  dynamic  programming  as  the  mathematical  optimization  technique. 
Although  dynamic  programming  is  very  powerful,  efficiency  deteriorates  rapidly  with  an  increase  in  the 
number  of  storage  tanks.  Since  OCOPS  has  a  water  distribution  simulation  routine  embedded  within  it, 
the  number  of  storage  tanks  that  can  be  included  in  the  optimization  is  limited  by  WADISO. 

Some  previous  methodologies  have  used  simplified  representations  or  approximations  of  the  distribu¬ 
tion  system  hydraulics.  In  these  approaches,  it  is  possible  for  error  to  be  introduced  into  the  optimization 
process.  In  OCOPS,  the  governing  equations  of  pressure  and  flow  distribution  are  solved  as  necessary 
and  no  simplifying  assumptions  are  made.  This  approach  allows  the  optimization  model  to  simulate 
extreme  operating  conditions. 

All  constraints  in  the  optimal  control  problem  formulated  in  this  study  are  assumed  to  be  binding 
and  cannot  be  violated.  Other  techniques  have  incorporated  penalty  functions  in  the  constraint  set.  Under 
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the  latter  approach,  constraints  can  be  violated.  However,  binding  constraints  prohibit  constraint  violation. 
Thus,  the  need  for  a  well  calibrated  model  of  the  distribution  system  is  readily  apparent  since  OCOPS 
initially  satisfies  all  constraints  and  then  attempts  to  optimize  the  system. 

OCOPS  is  very  flexible.  Proposed  changes  to  the  distribution  system  can  be  analyzed  by  simply 
modifying  the  computer  model  of  the  netrvork.  As  a  result,  the  impact  of  a  change  to  the  system  can  be 
quantified  and  used  to  help  decide  if  the  change  should  be  made. 

OCOPS  is  written  in  FORTRAN-77  and  can  be  run  on  any  IBM-compatible  PC  with  640K  of 
memory.  The  1  August  1988  simulation  took  nearly  8  hr  on  a  80386-based  computer.  Proposed  measures 
to  reduce  the  amount  of  computation  time  are  listed  in  Appendix  A. 

Most  computer  time  is  spent  computing  the  partial  derivatives  of  the  objective  function  and  of  the 
problem  constraints.  A  finite  difference  approximation  is  made,  requiring  a  number  of  WADISO  calls  for 
each  decision  variable  analyzed.  By  computing  the  derivatives  analytically,  computer  run  time  could  be 
significantly  reduced.  Furthermore,  analytically  computed  derivatives  eliminate  any  error  introduced  by 
the  finite  difference  approximation.  This,  in  turn,  may  cause  the  optimization  algorithm  to  move  closer 
to  a  global  optimum. 

The  OCOPS  model  should  be  reapplied  to  a  well  calibrated  system  to  fully  determine  the  model’s 
capabilities.  A  well  calibrated  model  of  a  water  distribution  system  can  be  created  using  high-quality  data 
on  system  tank  levels,  demands,  and  layout. 

Research  is  needed  on  analytically  computing  the  derivatives  of  the  objective  function  a r.d  constraints 
to  reduce  computation  times.  The  OCOPS  model  can  be  run  on  a  PC,  but  the  program  requires  a  large 
number  of  computations.  As  noted  above,  the  number  of  computations  could  be  reduced  substantially  by 
computing  the  derivatives  analytically.  Only  in  this  way  can  the  model  truly  approach  real-time  operation. 

OCOPS  should  be  modified  to  account  for  electricity  demand  charges.  At  present,  the  program 
considers  only  electricity  unit  charges.  For  some  water  utilities,  electricity  demand  charges  comprise 
the  major  fraction  of  electrical  costs.  These  additional  costs  would  be  incorporated  in  the  model  through 
additional  constraints. 

After  these  modifications  are  made  and  validated,  the  model  can  be  combined  with  a  demand 
forecasting  model  and  integrated  into  a  Supervisory  Control  and  Data  Acquisition  (SCADA)  system  to 
provide  automated  optimal  control  of  water  supply  pump  stations.  The  end  result  will  be  efficient 
operation,  reduced  electrical  costs,  and  wise  use  of  limited  natural  resources. 


METRIC  CONVERSION  TABLE 

1  hp  =  0.75  kW 
1  gal  =  3.785  L 
1  ft  =  0.305  m 
1  psi  =  6.895  kPa 
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APPENDIX  A: 


SOLUTION  METHODOLOGY 


Problem  Formulation 

Objective  Function 

The  optimal  control  problem  can  be  formulated  mathematically  as  a  nonlinear  optimization  problem 
subject  to  a  series  of  constraints.  In  the  approach  presented  in  this  report,  pump  operating  time  is  treated 
as  the  decision  variable.  The  operating  period,  usually  24  hr,  is  divided  into  a  series  of  T  time  intervals. 
For  a  water  distribution  system  with  I  pumps,  the  objective  function  of  the  optimal  control  problem  can 
be  written  as: 


Minimize 


I 

k 


[Eq  Al] 


where  Xti  is  the  duration  of  time  pump  i  operates  during  time  interval  t  FCX^)  is  the  total  cost  of 
pumping  for  pump  i  during  time  interval  t. 

Figure  Al  shows  a  typical  set  of  decision  variables  for  a  given  24  hr  period.  In  this  case,  the 
operating  period  is  divided  into  six  time  intervals  of  4  hr  each.  If  a  pump  operates  during  a  time  interval, 
it  is  assumed  to  begin  pumping  at  the  start  of  the  time  interval  and  operate  continuously  for  X  hr.  For 
instance,  during  the  hours  between  8:00  a.m.  and  12:00  p.m.,  pump  2  runs  for  2.5  hr  from  8:00  a.m.  to 
10:30  a.m.;  thus,  X3J  =  2.5. 

Problem  Constraints 

The  nonlinear  optimization  problem  is  subject  to  a  series  of  problem  constraints.  These  constraints 
can  be  categorized  as  (1)  explicit  bound  constraints  on  the  decision  variables,  (2)  implicit  bound 
constraints  on  system  parameters,  and  (3)  implicit  system  constraints. 

The  explicit  bound  constraints  set  lower  and  upper  limits  on  the  values  of  the  decision  variables.  For 
the  optimal  control  problem,  the  explicit  bound  constraints  can  be  expressed  as: 

0  <  X.i  <  At  [Eq  A2] 

where  At  is  equal  to  time  interval  t .  This  constraint  restricts  the  time  a  pump  can  operate  from  0  to  At 
hr  during  a  given  time  interval;  however,  a  pump  can  operate  during  any  or  all  time  intervals. 

The  implicit  bound  constraints  set  limits  on  system  parameters  such  as  pressure,  flow  rate,  system 
storage,  or  total  energy  consumption  for  a  given  time  interval.  For  any  lime  interval,  the  system  pressure 
at  any  node  may  be  bound  between  a  minimum  acceptable  pressure  and  a  maximum  allowable  pressure. 
The  implicit  pressure  constraint  is  expressed  as: 

Pmin*,  <  pk|  <  Pmaxk  (  [Eq  A3] 

where  Pk  t  is  the  pressure  at  node  k  during  time  interval  t.  High  pressures  (>  100  psi)  are  undesirable 
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Figure  Al.  Representation  of  the  decision  variables.  Shaded  areas  indicate  pump  is  running. 


since  they  have  a  tendency  to  waste  water  and  could  damage  residential  and  commercial  plumbing.  Low 
pressures  (<  40  psi)  do  not  allow  water  to  be  delivered  to  building  floors  above  four  stories.  Furthermore, 
many  health  codes  require  that  pressures  not  fall  below  20  psi.  Pressures  may,  however,  fluctuate  between 
Pmin  and  Pmax  during  any  time  interval. 

The  flow  rate  through  any  given  pipe  p  may  also  be  bound  between  an  acceptable  minimum  and 
allowable  maximum  during  a  given  time  interval.  This  constraint  can  be  written  as: 

Qmirip,!  <  Qp[  <  Qmaxp>t  [Eq  A4] 

where  Qp[  is  the  flow  rate  in  pipe  p  during  lime  interval  t.  Flow  rates  can  be  bounded  to  avoid  excessive 
head  loss  in  a  line.  For  this  study,  this  constraint  was  not  included  in  the  problem  formulation. 

Among  the  more  important  implicit  system  constraints  are  the  bounds  placed  on  tank  water  level. 
The  level  in  tank  m  can  be  bound  between  a  minimum  and  a  maximum  level.  The  tank  level  constraint 
can  be  expressed  as: 

Lminm,t  <  Lmil  <  Lmax™  t  [Eq  A5] 

where  Lm  t  is  the  water  level  in  tank  m  during  time  period  t.  In  most  cases,  the  minimum  allowable  tank 
level,  Lmir^j,  will  not  be  the  actual  bottom  of  the  tank.  Rather,  Lrninml  will  be  several  feet  above  the 
bottom  of  the  tank  to  allow  for  a  reserve  of  water  for  emergencies,  including  fire  fighting.  Tank  levels 
can  vary  between  Lmin  and  Lmax  during  any  time  interval.  It  is  the  energy  stored  in  the  tanks  that  allows 
pump  operation  to  be  scheduled  in  such  a  way  as  to  effect  least-cost  pumping. 

The  total  electrical  consumption  of  the  pump  station  during  any  time  interval  may  be  bounded 
between  a  lower  and  upper  value.  This  enables  control  of  kilowatt  consumption  so  that  electrical  demand 
charges  can  be  minimized.  The  constraint  for  electrical  consumption  can  be  written  as: 

Emir^  <  Et  <  Emaxt  [Eq  A6] 

where  E,  is  the  combined  electrical  consumption  of  all  pumps  subject  to  the  demand  charge  during  time 
interval  t.  Kilowatt  consumption  can  vary  between  Emin  and  Emax  during  a  time  interval. 

The  implicit  system  constraints  represent  the  equations  of  conservation  of  mass  and  energy  associated 
with  hydraulic  networks,  and  describe  the  state  of  the  system.  The  nodal  conservation  of  mass  constraint 
requires  that  the  sum  of  pipe  flows  into  or  out  of  a  node,  less  any  external  or  internal  demands,  be  equal 
to  zero.  A  conservation  of  mass  constraint  is  necessary  for  each  node  and  can  be  expressed  as: 


F(Q)  =  2*  Q  —  h  *0  k=  1 , 2 , 3 , . . .  K 

k  n“ 1  n  k 


|Eq  A7] 


where  Qn  is  the  flow  rate  in  pipe  n,  Dk  is  the  demand  at  node  k,  Nk  is  the  number  of  pipes  connected  at 
node  k,  and  K  is  the  number  of  nodes  in  the  system. 

Since  the  flow  in  a  pipe  is  a  function  of  the  head  difference  between  the  beginning  and  end  of  the 
pipe,  Equation  A7  can  be  rewritten  in  terms  of  head  difference.  This  essentially  combines  the 
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conservation  of  mass  and  energy  equations  into  a  single  equation.  The  advantage  of  combining  these 
equations  is  a  reduced  number  of  simultaneous  equations  that  must  be  solved.  When  using  the 
Hazen- Williams  head  loss  formula  under  steady-state  conditions,  the  combined  equation  can  be  expressed 
as: 
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where: 


k=l ,2 ,3, . .  .K  [Eq  A8] 
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Hb  =  Head  at  the  beginning  of  pipe  n,  ft 
He  =  Head  at  the  end  of  pipe  n,  ft 
Ln  =  Length  of  pipe  n,  ft 
Dn  =  Diameter  of  pipe  n,  ft 

Cn  =  Hazen- Williams  roughness  coefficient  for  pipe  n 
m  =  1.852. 


Optimization  Methodology 

To  solve  the  optimization  problem,  the  nonlinear  objective  or  cost  function  must  be  minimized  subject 
to  explicit  bound  constraints,  implicit  bound  constraints  on  system  parameters,  and  implicit  system  con¬ 
straints.  This  results  in  a  very  complex  problem.  The  complexity  of  the  problem  can  be  reduced 
substantially  by  linking  a  nonlinear  optimization  algorithm  with  a  water  distribution  simulation  routine. 
The  implicit  system  constraints  are  solved  by  the  simulation  routine  while  the  optimization  algorithm 
computes  the  decision  variables. 

That  was  the  approach  taken  in  the  optimization  methodology  presented  in  this  study.  The  nonlinear 
optimization  algorithm  used  was  GRG2.  The  simulation  routine  used  was  WADISO.  The  two  programs 
were  linked  together  using  the  GCOMP  subroutine  which  computes  the  values  of  the  objective  function 
and  the  constraints.  The  combined  algorithms,  along  with  a  main  program,  are  contained  in  the  computer 
program  OCOPS. 
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The  main  program  allows  data  entry  and  initiates  execution  of  the  nonlinear  optimization  algorithm. 
The  required  data  include  a  WADISO  system  file  and  an  optimization  data  file.  The  WADISO  system 
file  is  the  mathematical  model  of  the  distribution  system  and  contains  information  such  as  pipe  length, 
roughness,  diameter,  node  demand,  and  elevation.  Optimization  data  include  information  on  the  pumps 
and  tanks  in  the  system,  the  number  of  increments  in  the  operating  period,  the  initial  starting  point, 
electrical  costs,  and  GRG2  tolerances.  A  complete  listing  of  the  data  entry  portion  of  the  program  is 
provided  in  Appendix  B. 

The  GCOMP  routine,  which  evaluates  the  objective  function  and  problem  constraints,  is  required  by 
GRG2  to  determine  the  reduced  gradient,  step  size,  and  whether  a  solution  is  optimal.  This  routine  also 
updates  the  tank  levels.  For  the  optimal  control  problem,  the  objective  function  is  the  total  cost  of  pump¬ 
ing  for  a  24-hr  period. 

Accurately  computing  pumping  costs  over  a  24-hr  period  is  complicated  because  cost  is  a  function 
of  pump  head,  discharge,  and  efficiency.  To  further  complicate  matters,  pump  head  and  efficiency  vary 
with  discharge  and  pump  discharge  varies  with  system  demands  and  tank  water  levels.  Assuming 
electricity  costs  can  also  vary  over  time,  the  equation  required  to  compute  pumping  costs  for  one  pump 
can  be  represented  by: 

T 

Cost  -  A  t  Rt  ^  dt  [Eq  A9] 

JQ  t 


where: 

=  Pump  head  at  time  t,  ft 
Q,  =  Pump  discharge  at  time  t,  cu  ft/sec 
ft  =  Pump  efficiency  at  time  t 
Rt  =  Electrical  rate  at  time  t,  $/kWh 
T  =  Operating  period,  hr 
y  =  62.4  Ib/cu  ft. 


Pumping  costs  can  be  approximated  by  dividing  the  operating  period  T  into  several  time  intervals 
of  t  hr  each  and  finding  the  average  cost  over  each  time  interval.  That  is  the  approach  taken  by  OCOPS. 
At  the  beginning  of  each  time  interval,  the  implicit  system  constraints  are  satisfied  by  WADISO  and 
values  for  pump  head  and  discharge  are  passed  back  to  the  GCOMP  subroutine  where  the  cost  of  pumping 
is  computed.  Flows  in  lines  connecting  tanks  to  the  system  are  also  passed  to  the  GCOMP  routine  so  that 
tank  levels  can  be  updated.  At  the  end  of  the  time  interval,  the  system  is  rebalanced  based  on  the  updated 
tank  levels.  Again,  pump  head  and  discharge  values  are  passed  to  GCOMP  so  that  a  cost  of  pumping 
corresponding  to  the  new  system  conditions  at  the  end  of  the  time  interval  can  be  computed.  The  two 
pumping  costs  are  averaged  to  find  the  cost  of  pumping  for  the  time  interval. 
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During  any  time  interval,  pumps  can  be  taken  oilline  or  tanks  can  close  full  or  empty.  II  this  occurs, 
the  time  to  such  a  change  is  found,  the  implicit  system  constraints  arc  satisfied  immediately  before  the 
change,  values  for  pump  head  and  discharge  arc  passed  to  GCOMP,  tank  levels  are  updated,  and  an 
average  cost  for  the  subinterval  is  computed.  The  system  is  also  balanced  immediately  after  the  change. 
Critical  information  is  passed  to  GCOMP  where  it  is  used  to  compute  the  average  cost  for  the  subinterval 
corresponding  to  the  next  change  or  the  end  of  the  time  interval.  Tank  levels  are  also  updated  for  the  next 
subinterval. 


The  mathematical  form  of  the  objective  function  for  i  number  of  pumps  using  several  time  intervals 
can  be  expressed  as: 


Cost 


=  V  V  0^746! 

2,  L  550  t,i 


t=i  1=1 


,  i 


*  R  *  X„  , 
t  t  ,i 


where: 

Htj  =  Pump  head  of  pump  i  during  time  interval  t,  ft 
Qt ;  =  Pump  discharge  of  pump  i  during  time  interval  t,  cu  ft/sec 
C,  j  =  Efficiency  of  pump  i  during  time  interval  t 
Rt  =  Electrical  rate  during  time  interval  t,  $/kWh 
Xt j  =  Time  pump  i  runs  during  time  interval  t,  hr 


y  =  62.4  lb/sq  ft 
T  =  Operating  period,  hr 
I  =  Total  number  of  pumps  optimized. 


[Eq  A 10| 


In  addition  to  computing  the  objective  function,  the  GCOMP  subroutine  also  computes  the  values  of 
the  problem  constraints.  In  OCOPS,  the  problem  constraints  computed  by  GCOMP  are  pressure  at  critical 
nodes  and  the  tank  water  levels.  At  the  end  of  each  time  interval,  the  implicit  system  constraints  are 
solved  by  WADISO.  The  head  at  each  node  and  the  flow  in  each  line  connected  to  a  tank  are  passed  back 
to  the  GCOMP  routine.  This  information  is  used  to  find  the  pressure  at  critical  nodes  and  to  update  tank 
water  levels  so  that  values  for  the  problem  constraints  can  be  found.  A  flowchart  of  the  GCOMP  routine 
is  shown  in  Figure  A2  and  a  complete  listing  of  the  subroutine  is  provided  in  Appendix  C. 
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GRG2  Solution  Algorithm 


GRG2  is  an  optimization  algorithm  that  solves  nonlinear  problems  using  a  generalized  reduced 
gradient  technique.  Problems  are  of  the  general  form: 

Minimize  F(X)  [Eq  All] 

Subject  to: 

LBj  <  Gj(X)  <  UBj  i=l,2,3...m  [Eq  A12] 

LBj  <  Xj  <  UBj  j=l,2,3...n  [Eq  A13] 

X  is  a  one-dimensional  vector  with  n  elements  called  decision  variables.  F(X)  is  the  objective 
function  that  is  minimized  subject  to  the  set  of  problem  constraints,  Gj(X).  There  are  m  linear  or 
nonlinear  constraints  to  the  problem,  each  with  a  lower  bound,  LBj,  and  an  upper  bound,  UBj.  It  is 
assumed  that  there  are  more  decision  variables  than  problem  constraints;  otherwise,  the  problem  would 
be  infeasible  or  have  a  unique  solution.  There  are  also  lower  and  upper  bounds  on  the  decision  variables. 

The  underlying  theory  of  GRG2  is  to  find  a  solution  for  m  of  the  variables,  or  basic  variables,  in 
terms  of  the  remaining  n-m  variables,  or  nonbasic  variables.  This  is  accomplished  through  rewriting  the 
constraint  set  as  equality  constraints  by  adding  slack  variables.  The  new  constraints  written  in  terms  of 
basic  and  nonbasic  variables  are  expressed  as: 

gi(y(x),x)  =  0  i=l,2,3...m  [Eq  A14] 


where: 

y(x)  =  Vector  of  basic  variables 
x  =  Vector  of  nonbasic  variables. 

The  nonlinear  problem  can  now  be  transformed  into  a  reduced  problem  expressed  as: 

MIN  F(X)  =  f(y(x),x)  [Eq  A 15] 

subject  only  to  lower  and  upper  bounds  on  the  nonbasic  variables.  GRG2  solves  the  original  nonlinear 
problem  by  solving  a  series  of  reduced  unconstrained  minimization  problems. 

GRG2  can  optimize  problems  if  an  infeasible  starting  point  is  provided.  An  infeasible  starting  point 
is  one  where  at  least  one  problem  constraint  is  violated.  In  this  case,  the  algorithm  enters  into  a  Phase  1 
optimization  where  the  objective  function  becomes  the  sum  of  the  constraint  violations.  The  Phase  1 
objective  function  is  minimized  until  all  constraints  arc  within  their  bounds.  At  this  point,  the  Phase  2 
optimization  begins. 

Once  the  reduced  problem  has  been  formulated  and  a  feasible  starting  point  obtained,  GRG2  uses  the 
following  steps  to  minimize  the  objective  function. 
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1.  The  basis  inverse  and  reduced  gradient  are  computed  and  used  to  test  the  optimality  of  the  feasible 
starting  point.  Optimality  is  checked  via  the  Kuhn-Tucker  conditions  or  if  the  fractional  change  in  the 
objective  function  is  less  than  a  specified  tolerance  for  a  specified  number  of  iterations. 

2.  If  the  starting  point  is  not  optimal,  a  new  decision  vector  (xnew)  is  generated  using  the  following 
recursive  equation: 

Xnew1  =  Xoldi  +  a  di  1  =  l.2,3...n  [Eq  A16] 

where  dj  is  the  search  direction  and  a  is  the  step  size. 

3.  The  search  direction  is  computed  by: 

di  =  -H.  \7F(X  )  i  =  1  ,2,3. .  .n  [Eq  A17] 

where  -Hj  is  an  n  by  n  symmetric  positive  definite  matrix  and  VF(X,)  is  the  reduced  gradient. 

4.  The  magnitude  of  a  is  determined  by  substituting  Equation  A16  into  Equation  All  and 
determining  the  value  of  ot  that  minimizes  the  resulting  one-dimensional  objective  function  shown  below: 

MIN  F  (Xi  +  a  d.)  [Eq  A18] 

5.  The  optimal  value  of  the  step  size,  a,  is  determined  using  a  one-dimensional  line  search.  During 
the  line  search,  constraints  on  the  original  nonlinear  problem  are  checked  for  violation.  If  a  constraint  is 
violated,  the  step  size  will  be  adjusted  to  keep  the  search  within  the  bounds  of  the  feasible  region.  Once 
the  optimal  value  for  a  is  found,  the  new  decision  vector  is  determined  using  Equation  A 16  and  the 
process  is  repeated. 

Note:  if  a  variable  crosses  an  explicit  bound,  that  variable  is  set  to  the  value  of  the  bound  and  a 
change  of  basis  occurs.  The  change  of  basis  results  in  the  bounded  variable  leaving  the  basis  and  a  new 
variable  not  at  a  bound  entering  the  basis.  At  this  point,  optimality  is  checked.  If  the  new  point  is  not 
optimal,  the  process  is  repeated. 

WADISO  Solution  Algorithm 

WADISO20  is  a  computer  simulation  model  used  to  determine  the  pressure  and  flow  distribution 
within  a  water  distribution  system.  Flow  and  pressure  within  a  pipe  network  can  be  determined  by  either 
of  two  methods:  (1)  the  loop  method  or  (2)  the  node  method.  The  loop  method  initially  solves  for  the 
flow  in  the  pipes  and  backs  out  the  nodal  heads.  In  contrast,  the  node  method  solves  for  the  heads  at  the 
nodes  and  then  computes  pipe  flow.  WADISO  uses  the  node  method  to  determine  the  pressure  and  flow 
distribution  in  a  network.  Both  methods  provide  accurate  results;  however,  the  node  method  offers  several 
advantages  over  the  loop  method  in  the  optimal  control  application  presented  here. 


“J.  Gesslcr  and  T.  M.  Walski. 
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The  number  of  equations  that  must  be  solved  to  obtain  the  pressure  and  flow  distribution  is  generally 
less  with  the  node  method.  This  is  especially  true  for  large  systems  since  the  number  of  equations  is  equal 
to  the  number  of  nodes.  In  the  loop  method,  the  number  of  equations  is  equal  to  the  number  of  pipes  in 
the  system.  Therefore,  computer  storage  and  execution  lime  can  be  reduced  by  using  the  node  equations. 
This  choice  of  method  becomes  very  important  in  optimizing  large  pipe  networks  because  the  equations 
describing  the  water  distribution  system  must  be  solved  many  times. 

The  loop  method  requires  that  the  system  geometry  be  determined  so  that  an  energy  path  between 
two  known  head  sources  can  be  found.  Analysis  of  the  system  geometry  is  not  necessary  with  the  node 
method.  Also,  due  to  the  nature  of  the  node  equations,  pressure  control  devices  such  as  pumps,  pressure 
reducing  valves,  and  check  valves  can  be  included  in  the  computer  model  without  any  increase  in 
computer  time. 

The  node  method  is  based  on  the  equation  of  continuity  which  states  that  all  mass  arriving  at  a  node 
(Qin)  must  equal  all  mass  leaving  a  node  (Qout  +  D).  Fluid  can  reach  a  node  as  flow  through  other 
connecting  pipes  or  as  a  system  demand  placed  at  the  node.  Continuity  at  a  node  is  expressed  as: 
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[Eq  A19] 


The  flow  rate  in  a  pipe  is  a  function  of  the  head  difference  between  the  beginning  and  end  of  a  pipe 
as  shown  by  the  Hazen-Williams  formula- 
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[Eq  A20] 


As  a  result,  the  equation  of  continuity  can  be  rewritten  in  terms  of  the  head  drop  across  a  pipe  and  the 
characteristic  pipe  coefficient  as  shown  below: 
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[Eq  A21] 


WADISO  also  can  incorporate  pumps  in  the  system  equations.  A  pump  is  considered  to  be  a  pipe 
with  a  head  flow  relationship  equal  to  the  pump  curve.  The  head-flow  characteristics  of  a  pump  can  be 
closely  approximated  by  a  parabola  of  the  form: 

Hp  =  A  +  BQ  +  CQ2  [Eq  A22] 

Using  the  quadratic  rule  to  solve  a  second-order  polynomial  and  further  assuming  that  the  linear  and 
squared  coefficients  are  negative,  the  discharge  from  a  pump  can  be  expressed  as  a  function  of  head  as 
shown  in  Equation  A23.  It  is  necessary  that  the  linear  term  of  the  pump  characteristic  curve  be  negative 
since  a  positive  value  would  imply  a  curve  with  a  maximum  head  at  a  positive  discharge.  A  negative 
squared  term  ensures  that  the  curve  is  concave.  Either  of  these  conditions  could  cause  convergence 
problems  while  solving  the  system  of  equations. 
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If  one  equation  expressing  mass  continuity  is  written  for  each  node  in  the  system,  then  a  K  by  K  (K 
equals  the  number  of  nodes  in  the  system)  coefficient  matrix  of  node  equations  can  be  established. 
Fortunately,  the  coefficient  matrix  is  symmetrical  and  sparse  so  that  only  half  of  the  matrix  is  used  to 
solve  the  equations,  thereby  reducing  computer  time.  WADISO  uses  a  Gaussian  elimination  technique 
to  solve  the  equations  simultaneously.  Determining  the  pressure  and  flow  distribution  is  an  iterative 
process  with  the  flow  rates  of  the  previous  iteration  updating  the  coefficient  matrix  of  the  current  iteration. 


General  Considerations 

The  optimal  control  problem  is  very  large  and  complex.  The  size  of  the  problem  and  the  computer 
time  required  to  solve  it  can  be  reduced  by  adhering  to  several  simple  guidelines. 

The  number  of  decision  variables  is  equal  to  the  product  of  the  number  of  pumps  to  be  optimized 
and  the  number  of  time  intervals  in  the  24-hr  time  period.  The  computer  time  necessary  to  solve  the 
optimal  control  problem  can  be  reduced  significantly  by  reducing  the  number  of  decision  variables.  These 
variables  can  be  reduced  by  including  only  those  pumps  run  on  a  daily  or  weekly  basis.  In  other  words, 
any  pumps  that  are  usually  on  stand-by  or  used  only  in  case  of  emergency  are  not  included  in  the  analysis. 
Also,  pumps  that  use  an  alternative  fuel,  such  as  diesel  fuel,  should  not  be  included  in  the  analysis  since 
the  methodology  assumes  all  pumps  are  driven  by  electricity. 

Another  way  to  reduce  the  size  and  complexity  of  the  problem  is  to  skeletalize  the  water  distribution 
system.  When  using  a  computer  model  of  a  water  distribution  system,  it  is  not  necessary  to  include  all 
pipes  in  the  system.  Rather,  only  the  major  pipelines  are  included.  Results  from  the  skeletal  model  are 
still  accurate  because  those  lines  not  included  do  not  carry  a  large  amount  of  flow.21 

WADISO  uses  the  nodal  method  to  solve  the  equations  of  continuity  and  energy  associated  with 
water  distribution  systems.  Therefore,  by  reducing  the  number  of  pipes  in  the  mathematical  model  of  the 
system,  the  number  of  nodes  can  be  lessened  and  computation  time  significantly  reduced. 

In  most  distribution  systems,  pumps  are  placed  online  at  any  time  during  the  day.  In  other  words, 
pumps  operate  over  a  continuous  timeframe.  In  the  optimal  control  methodology  presented  here,  pumps 
operate  over  a  discrete  timeframe  since  it  is  assumed  that  pumps  are  placed  online  at  the  beginning  of  a 
time  interval.  The  smaller  the  time  interval  used,  the  closer  the  approximation  to  continuous  operation. 

Ideally,  a  very  small  time  interval  would  be  used  since  it  affords  more  flexibility  in  pump  operation. 
Unfortunately,  the  size  of  the  time  interval  is  inversely  proportional  to  the  computer  run  time  and  required 
computer  storage.  Computer  run  time  can  be  reduced  substantially  by  using  a  larger  time  step  equal  to 
3  or  4  hr  as  opposed  to  0.5  or  1  hr.  However,  increased  savings  in  pump  cost  may  be  realized  by  using 
a  smaller  time  interval. 
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GRG2  requires  an  initial  guess  to  start  the  optimization  process.  The  guess  can  be  either  feasible  or 
infeasible.  A  feasible  initial  guess  is  one  that  does  not  violate  any  problem  constraints  and  an  infeasible 
guess  is  one  that  violates  at  least  one  problem  constraint.  GRG2  can  optimize  pump  operation  even  if  an 
infeasible  starting  point  is  provided;  however,  the  time  required  to  solve  the  problem  is  increased 
substantially  by  such  a  point. 

Several  simple  checks  can  be  made  to  determine  if  an  infeasible  starting  point  is  used.  If  a  pump 
operating  time  is  less  than  zero  or  greater  than  the  specified  time  interval,  GRG2  will  change  the  value 
of  the  decision  variable  to  its  lower  or  upper  bound,  respectively.  Although  this  type  of  infeasible  starting 
point  has  very  little  effect  on  run  time,  it  is  good  practice  to  provide  an  initial  time  of  operation  between 
zero  and  the  specified  time  interval. 

Node  pressures  have  been  formulated  as  problem  constraints.  The  problem  size  can  be  reduced  by 
specifying  only  critical  nodes  as  pressure  constraints.  In  other  words,  only  examine  those  nodes  most 
likely  to  have  pressures  that  violate  an  upper  or  lower  boundary.  Nodes  at  high  or  low  elevations  or  those 
at  the  extreme  end  of  the  distribution  system  can  be  considered  critical  nodes. 

If,  for  a  given  initial  guess,  the  pressure  constraint  is  violated,  then  a  new  initial  guess  that  does  not 
violate  the  constraint  can  be  used.  The  new  guess  should  reflect  how  and  when  the  constraint  was 
violated,  and  be  based  on  knowledge  of  the  system.  For  example,  if  the  pressure  at  a  node  is  below  the 
acceptable  minimum  during  time  step  3,  then  a  pump  can  be  placed  online  during  this  or  a  previous  time 
step. 


Tank  water  level  at  the  end  of  the  24-hr  period  has  also  been  formulated  as  a  problem  constraint. 
This  was  done  because  it  is  desirable  to  have  tank  levels  at  a  predetermined  elevation  at  the  end  of  the 
operating  period.  The  problem  has  been  formulated  so  that  the  specified  elevation  is  equal  to  the  lower 
bound  of  the  constraint.  The  nature  of  the  optimal  control  problem  forces  the  tank  levels  toward  their 
lower  bound.  As  a  result,  the  specified  ending  elevation  is  met. 

If  an  initial  guess  provided  is  such  that  at  least  one  lank  level  at  the  end  of  the  operating  period  is 
less  than  the  specified  elevation,  a  problem  constraint  is  violated.  By  turning  on  pumps  or  running  pumps 
longer  during  the  last  time  interval,  the  tank  constraint  may  not  be  violated.  However,  GRG2  can 
optimize  problems  with  an  infeasible  starting  point. 
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APPENDIX  B: 


MAIN  DATA  ENTRY  ROUTINE 


The  main  data  entry  routine  opens  and  reads  information  from  two  data  files--one  describing  the 
water  distribution  system  and  another  containing  information  on  pump  characteristics,  electrical  costs,  etc. 
The  file  that  describes  the  water  distribution  system  is  actually  a  WADISO  data  file  created  using  the 
WAD1SO  computer  program.  A  WADISO  system  file  is  the  mathematical  representation  or  computer 
model  of  the  water  distribution  system.  Information  contained  in  this  file  includes  pipe  length,  diameter, 
roughness,  node  demand,  and  elevation.  The  pump  head-discharge  relationship  is  also  contained  in  this 
file. 


The  second  file  opened  in  the  main  data  entry  routine  defines  the  number  of  pumps  and  tanks  in  the 
system  and  the  time  interval  to  be  used  for  the  analysis.  The  relationship  between  pump  efficiency  and 
pump  discharge  is  also  contained  in  this  file.  Other  information  includes  the  initial  starting  point  for  the 
nonlinear  optimization  routine,  the  cost  of  electrical  energy,  tank  geometry  and  connecting  pipes,  the 
global  demand  factor  for  each  time  increment,  and  GRG2  tolerances. 

The  main  data  entry  routine  also  defines  variables  and  initializes  several  arrays  used  by  GRG2. 
Variables  defined  by  this  routine  include  the  number  of  decision  variables  and  the  number  of  problem 
constraints.  The  arrays  include  the  initial  guess  of  the  decision  variables,  the  upper  and  lower  bounds  on 
the  decision  variable,  and  the  upper  and  lower  bounds  on  the  problem  constraints.  The  nonlinear 
optimization  algorithm  is  called  from  the  main  data  entry  routine. 

The  specific  format  for  the  WADISO  system  file  can  be  found  in  the  WADISO  user’s  guide.22  The 
WADISO  program  and  documentation  can  be  obtained  from  the  author.  The  formal  of  the  system 
optimization  data  file  is  described  in  Table  B1  and  a  sample  input  file  is  shown  in  Figure  Bl.  The  sample 
data  file  assumes  a  4-hr  time  interval  with  four  pumps  in  the  distribution  system  available  for  optimization. 
The  source  code  for  the  main  data  entry  routine  is  provided  in  Figure  B2. 
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Table  B1 


Format  of  System  Optimization  Data  File 


Variable 

Card  Number 

Format 

Column 

Description 

Name 

1 

15 

1-5 

Number  of  Pumps  Available  for 

Optimization 

NPUMP 

15 

6-10 

Time  Step  Used  (Hr) 

NSTEP 

15 

11-15 

Number  of  Tanks  in  the  System 

NTANK 

15 

16-20 

Print  Flag 

IBUG 

Pump  Efficiency  Card 

One  Card  for  Each  Pump  Optimized:  1=1,  NPUMP 


2 

15 

1-5 

- - - - r - r-.? - >  •  . .  . 

Link  Number  of  Pump 

LINK(I) 

F10.2 

6-15 

First  Discharge  Point 

X0 

FI  0.2 

16-25 

Efficiency  Corresponding  to  First 

Discharge  Point 

Y0 

FI  0.2 

26-35 

Second  Discharge  Point 

XI 

FI  0.2 

3645 

Efficiency  Corresponding  to  Second 

Discharge  Point 

Y1 

F10.2 

46-55 

Third  Discharge  Point 

X2 

FI  0.2 

56-65 

Efficiency  Corresponding  to  Third 

Discharge  Point 

Y2 

Initial  Starting  Point  Card 

One  Card  for  Each  Decision  Variable;  1=1,  NVAR 

3 

FI  0.2 

1-10 

Initial  Time  of  Operation  for  Pump  I 

During  Time  Interval  (Hr) 

GUESS 

FI  0.2 

11-20 

Electrical  Rate  in  Effect  for  This  Pump 

During  This  Time  Interval  (Cents/kWh) 

ERATE(I) 

Node  Constraint  Data  Card 

4 

15 

1-5 

Number  of  Critical  Nodes  in  the  System 

NNODES 

FI  0.2 

6-15 

Minimum  Acceptable  Pressure  (psi) 

PM  IN 

F10.2 

16-25 

Maximum  Allowable  Pressure  (psi) 

PMAX 

5 

1615 

1-5 

Node  Number  of  Critical  Node 

NNNODEd) 

1=1,  NNODES 
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Table  B1  (Cont’d) 


Card  Number 

Format 

Column 

Description 

Variable 

Name 

Tank  Data  Card 

One  Card  for  Each  Tank;  1=1,  NTANK 

6 

15 

1-5 

Node  Number  of  Tank 

NNNODEd) 

15 

6-10 

Link  Number  of  Tank  Connecting  Line 

NCON(I) 

F10.2 

11-20 

Top  Elevation  of  Tank  (ft) 

TOPELd) 

FI  0.2 

21-30 

Bottom  Elevation  of  Tank  (ft) 

BOTEL(I) 

FI  0.2 

31-40 

Tank  Diameter  (ft) 

DIAM 

F10.2 

41-50 

Tank  Level  Desired  at  the  End  of  the 

24 -Hr  Period  (ft) 

ENDELfl) 

FI  0.2 

51-60 

Minimum  Allowable  Tank  Level  for 

Fire  Storage  (ft) 

STOMIN(I) 

System  Demand  Data  Card 


7 

8F10.2 

1-10 

Global  Demand  Factor  for  This  Time 

Interval  1=1.  (24/NSTEP) 

RATIOd) 

GRG2  Data  Card  No.  1 

8 

A90 

1-90 

Title  for  This  Simulation 

TITLE 

GRG2  Data  Card  No.  2 

9 

F15.14 

1-15 

A  Constraint  Is  Assumed  To  Be  Binding 

If  It  Is  Within  This  Value  (FPNEWT) 
of  One  of  Its  Bounds 

FPNEWT 

F15.14 

16o0 

If  I.  Is  Desired  ^  Run  the  Problem 

With  FPNEWT  Initially  Set  Fairly  Large 
and  Then  Tightened  at  the  End  of  the 
Optimization,  This  Is  Accomplished 
by  Assigning  FP1N1T  the  Initial  Toler¬ 
ance  and  FPNEWT  the  Final  One 

FPINT 

FI  5.14 

31-45 

Minimum  Fractional  Change  in  Objective 

Allowed  To  Avoid  Program  Termination 

FPSTOP 

F10.2 

46-55 

Step  Size  for  Finite  Difference  Approx¬ 
imation  of  Partial  Derivatives  of  Con¬ 
straints  and  Objective  Function 

PSTEP 
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Tabic  B1  (Cont’d) 


Card  Number 

Format 

Column 

Description 

Variable 

Name 

GRG2  Data  Card  No.  3 

10 

15 

1-5 

If  Fractional  Charge  in  Objective 

Function  Is  Less  Than  FPSTOP  for 

NNSTOP  Iterations,  Then  Program 

Terminates 

NNSTOP 

15 

6-10 

If  Newton  Routine  Requires  More  Than 

IITL1M  Iterations  Before  Converging, 

Then  Iterations  Are  Stopped  and  Cor¬ 
rective  Action  Is  Taken 

IITLIM 

15 

11-15 

Print  Flag 

IIPR 

15 

16-20 

If  IIPN#  is  Greater  Than  Zero,  Then 

IPR  Will  Be  Set  to  #  After  IIPN# 

Iterations 

IIPN# 

15 

21-25 

Method  for  Initial  Estimates  of 

Variables  for  Each  One-Dimensional  Search 

IIQUADBaac 

15 

26-30 

Method  for  Obtaining  Partial  Deriva¬ 
tive  of  Constraints  and  Objective 

Function 

0  -  Forward  Difference  Approximation 

1  -  Central  Difference  Approximation 

2  -  User-Supplied  Subroutine 

LDERIV 

15 

31-35 

Method  for  Determining  Conjugate 

Gradient 

MMODCG 
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4 

4  6 

1 

7 

800. 

69. 

1300. 

83. 

1800. 

74. 

8 

1600. 

71. 

2300. 

82. 

3000. 

75. 

9 

2000. 

72. 

3500. 

85. 

4700. 

70. 

11 

4000. 

78. 

5500. 

86. 

7500. 

80. 

0.00 

4.  50 

1 

7 

1 

4.00 

4.  50 

8 

*“> 

4m 

4.00 

4.  50 

9 

3 

0.  00 

4.  50 

11 

4 

1.60 

4.  50 

2 

7 

5 

4.00 

4.  50 

8 

6 

4.  00 

4.  50 

9 

7 

0.38 

4.  50 

11 

8 

4.00 

4.  50 

3 

7 

9 

4.  00 

4.  50 

8 

10 

4.00 

4.50 

9 

11 

1.  70 

4.  50 

11 

12 

4.00 

4.  50 

4 

7 

13 

4.  00 

4.  50 

8 

14 

4.00 

4.  50 

9 

15 

1.73 

4.50 

11 

16 

4.00 

4.  50 

5 

7 

17 

4.  00 

4.  50 

8 

18 

4.  00 

4.  50 

9 

19 

1.  43 

4.  50 

11 

20 

4.00 

4.  50 

6 

7 

21 

4.  00 

4.  50 

8 

22 

4.00 

4.50 

9 

23 

0.  24 

4.50 

11 

24 

3  30.0  130.0 


24 

42 

47 

25 

27 

37.34 

0.0 

47.74 

37.340 

9.  335 

32 

41 

34.75 

0.0 

69.99 

34.710 

8.688 

39 

49 

39.  17 

0.0 

46.61 

39. 170 

9.  790 

45 

57 

35.  00 

0.0 

85.  41 

35. 000 

8.  750 

50 

63 

39.  17 

0.0 

46.  61 

35. 000 

9.790 

57 

74 

35.  00 

0.0 

69.74 

21. 400 

8.668 

1.93581 

2. 56606 

3.  10616 

2.75870 

2.  88290 

1.98056 

FORT 

HOOD 

PUMP  OPTIMIZATION 

DATA  1  AUG  88 

0. 000000001 

.0001 

0. 

001 

3 

15 

3  0 

0  0 

1 

Figure  Bl.  Sample  input  file  for  the  main  data  entry  routine. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


»*««*»*«*»«***•*««*»**«*»**«*«****««***«»********«******** 

US  ARMY  COPRS  OF  ENGINEERS  PUMP  OPTIMIZATION  PROGRAM 
FOR  USE  TO  DETERMINE  LEAST  COST  PUMP  OPERATION 
GIVEN  A  VARIABLE  ELECTRIC  RATE  SCHEDULE,  SYSTEM 
DEMANDS,  AND  PUMP  CHARACTERISTICS. 

WRITTEN  BY:  DONALD  V.  CHASE 
MARCH/APRIL  1988 
US  ARMY  CORPS  OF  ENGINEERS 
WATERWAYS  EXPERIMENT  STATION 
CEWES-EE-R 
P.O.  BOX  631 
VICKSBURG,  MS  39180 
601/634-3931 

**#**##»#»#*#»»##»*»****###*#**»*«»#«**»**»#**#**»**#*»*** 

***####*#*«**»**##**#***#***#*#*#*##***#*#«»»##*»»»*#»**#* 


THIS  PROGRAM  WILL  CALL  GRG2,  A  NONLINEAR  OPTIMIZATION  ROUTINE, 
AND  OPTIMIZE  PUMP  OPERATION.  DECISION  VARIABLES  WILL  BE  TIME 
OF  OPERATION  OF  EACH  PUMP,  AND  CONSTRAINTS  WILL  BE  WATER  LEVEL 
IN  THE  SYSTEM  TANKS. 


IMPLICIT  REAL*8( A-H, 0-Z) 
IMPLICIT  INTEGER*4(I-N) 


CHARACTER  FILEIN*14,  FILE0UT*14,  TTITLE»90,  ST*90,  DATFLE»14 
CHARACTER  FILEDT*14, FILEAT*14 


DIMENSION  FCNSQ50),  RMULTS<150),  NONBASC 150),  REDGR(  150),  INBIND(150> 
DIMENSION  RAMCON( 150),  RAMVARt 150),  DEFAULt 19),  XX ( 150) 

DIMENSION  STOM I N ( 20 ) , Z ( 20000 ) 

INTEGER  PNL,  0 


PARAMETER  < PNL=200,  LNL=150,  MNL=20, IA=2000> 
PARAMETER  <NDC=150, NTP=20, NPP=10) 


COMMON 

COMMON 

COMMON 

1LFLAG( 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

ID 


/BLK1/  NVAR,  NTANK,  NSTEP,  NNOBJ,  Cl,  NN1,  NPUMP,  CDEM( 100) 

/BLK2/  ERATE ( NDC ) ,  AREA ( NTP ) ,  COEFF( NPP, 3 ) , LINK ( NPP ) 

/BLK2A/  CTANK(NTP), NNODE( NTP),  NCON( NTP),  RATIO (24), 

NPP), NFLAG(NTP) 

/BLK3/  XO,  YO,  XI,  Yl,  X2, Y2 
/BLK4/  C0E1,  C0E2,  C0E3 

/BLK5/  BLVAR(NDC),  BUVAR(NDC),  BLCON < 2«NDC ) , BUC0N(2»NDC) 
/TANK/  TOPEL(NTP), BOTEL< NTP ) , ENDEL( NTP ),  LTANK 
/NODES/  0(15), HE< PNL),  DO(PNL),  EL(PNL),  S( PNL),  Gl( PNL) 

/PIPES/  A< IA),  CP (PNL),  DI (PNL),  XL(PNL),  HW(PNL) 

/TOPOL/  IBEtPNL';  IEN(PNL),  IPI ( PNL) ,  INO(PNL) ,  IBI (PNL) ,  IEI ( PN 


Figure  B2.  Source  code  for  the  main  data  entry  routine. 
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non  on 


COMMON  /PRINT/  IPM,  IPP,  IPE,  NDFLG,  IBUG 
COMMON  /JOB/  JOB 

COMMON  / ACCU/  PRAC, FLAC, HWMA, ICL 
COMMON  /FILE/  IIN, IOUT 
COMMON  /BLK6/  PSTEP 

COMMON  /BLK8/  NNODES, PMIN, PMAX, NNNODE( 20), NN2, NN3 
LOGICAL  INPRNT, OTPRNT 


C 

C 

C  FILE  MANIPULATION  BLOCK 

C 

C 

C  WRITEU,  1001) 

C  READ(*,1002)FILEIN 

C  WRITE( *, 1003) 

C  READ(«, 1002)FILE0UT 

C 

C 

IIN  =  5 
I0UT=7 


OPEN ( IIN,  ERR=3060,  FILE= ' SYSTEM.  DAT' ,  STATUS= ' OLD ' , RECL=63 > 
OPEN( IOUT, FILE= ' OPTIMAL.  OUT' ,  STATUS= ' UNKNOWN ' ) 

REWIND  1 
REWIND  IIN 
REWIND  IOUT 

FILEIN= ' SYSTEM. DAT' 


GOTO  10 

3060  WRITE( IOUT, 1004) 

STOP 

10  CALL  SIMULATE(FILEIN) 

NDFLG=1 

C 

C  ***»*«*****»*###»»***»»###»*»»»##*»*#*****#****##»#*#*»*# 

C  HlimHIHXOIMXliHHHtIKOMHilOHKIHHIOmO 

C 

C  ENTER  INFORMATION  ON  DECISION  VARIABLES 

C 

C  DECISION  VARIABLES  FOR  PUMP  OPTIMIZATION  IS  TIME. 

C  THE  NUMBER  OF  DECISION  VARIABLES  IS  THE  PRODUCT  OF  THE 

C  NUMBER  OF  PUMPS  WHICH  CAN  BE  OPERATED  AND  THE  NUMBER  OF 

C  TIME  STEPS  SPECIFIED.  FOR  EXAMPLE,  IF  FIVE  (5)  PUMPS 

C  ARE  AVAILALBE  TO  THE  SYSTEM  AND  THE  TIME  STEP  SPECIFIED 

C  IS  FOUR  (4)  HOURS,  THE  NUMBER  OF  DECISION  VARIABLES  IS 

Figure  B2.  (Cont’d). 
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n  o 


C  30.  30  =  5  PUMPS  *  (24  HOURS/  4  HOURS) 

C 

C 

C  **#*«#»**»**#»**#***»***#*##**#»******#«»»#*»**####»»#**#» 

C  »*»**###»**#*****#*#***##»**#*##*«#*#*##»##»*#»»#*»#####** 

C 

C  INITIALIZE  VALUES  FOR  GRG  ROUTINE 

C 
C 

C  INPRNT : 

C 
C 
C 

C  OTPRNT: 

C 
C 
C 

C  NCORE: 

C  NVAR : 

C  NFUN : 

C  MAXBAS : 

C  MAXHES: 

C  NNOBJ : 

C  TTITLE : 

C  BLVAR: 

C  BUVAR : 

C  BLCON : 

C  BUCCN : 

C  DEFAUL: 

C 

C  #**»»##*#»*#»*#«**#*»#»#»*#»*##**##*#»#«###**»*#*#*#*##»*# 

c  ***#*##*#***#**#*#****#*»#*#***«#****»#*»*#****#******#»«# 

L00P=1 

C  WRITE!*, 1035) 

C  READ( *,  1002)FILEDT 

C  WRITE!*, 1036) 

C  READ!*, 1002)FILEAT 

OPEN ( UNIT=1,  FILE= 'OPTIMAL.  DAT',  STATUS* 'OLD' ) 

OPEN ( UNIT=8,  FILE* ' OUT ' ,  STATUS= ' UNKNOWN ' ) 


C  NPUMP  IS  THE  NUMBER  OF  PUMPS  AVAIALBE  FOR  OPERATION 

C  NSTEP  IS  THE  TIME  STEP  SPECIFIED  IN  HOURS.  A  24  HOUR 

C  ANALYSIS  IS  ASSUMED. 

C  LPFLG  IS  A  LOOP  FLAG 

NTANK  IS  THE  NUMBER  OF  TANKS  IN  THE  SYSTEM  CONTRIBUTING  TO 
HEAD  SEEN  BY  THE  PUMP  STATION 


LOGICAL,  ENABLES  PRINT  OPTION 
TRUE:  PRINT  INPUT  DATA 
FALSE:  DO  NOT  PRINT  INPUT  DATA 

LOGICAL,  ENABLES  PRINT  OPTION 
TRUE:  PRINT  FINAL  RESULTS 
FALSE:  DO  NOT  PRINT  FINAL  RESULTS 

DIMENSION  OF  THE  Z  ARRAY 

NUMBER  OF  DECISION  VARIABLES 

NUMBER  OF  FUNCTIONS  INCLUDING  OBJECTIVE 

UPPER  LIMIT  ON  THE  NUMBER  OF  BINDING  CONSTRAINTS 

MAXIMUM  ALLOWABLE  SIZE  OF  THE  HESSION  MATRIX 

INDEX  OF  THE  OJECTIVE  FUNCTION 

ALPHANUMERIC,  TITLE  OF  PROBLEM 

ARRAY  OF  LOWER  BOUND  OF  VARIABLES 

ARRAY  OF  UPPER  BOUND  OF  VARIABLES 

ARRAY  OF  LOWER  BOUND  OF  CONSTRAINTS 

ARRAY  OF  UPPER  BOUND  OF  CONSTRAINTS 

ARRAY  OF  DEFAULT  VALUES 
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C  PRAC  IS  THE  PRESSURE  TOLERANCE  USED  BY  WADISO 
C  FLAC  IS  THE  FLOWRATE  TOLERANCE  USED  BY  WADISO 
C  ICL  IS  THE  ITERATION  TOLERANCE  USED  BY  WADISO 
C  I3UG  IS  A  PRINT  FLAG 


READ< 1, 1005 ) NPUMP, NSTEP,  LPFLG, NTANK, PRAC, FLAC, ICL, TBUG 

C3=448. 831 

FLAC=FLAC/C3 

HWMA=100. 


NVAR  IS  THE  TOTAL  NUMBER  OF  DECISION  VARIABLES 

INTIALIZE  CONSTANTS 

Cl  =  PUMP  COEFFICIENT 

NN1  =  NUMBER  OF  TIME  STEPS 

Cl=<62.4»0.746)/550. 

NN1=24/NSTEP 

NVAR=NPUMP*NN1 

ECHO  INPUT 

WRITE ( IOUT,  1006 ) NPUMP,  NSTEP, NVAR, NTANK 


INTIALIZE  ARRAYS 
DO  100  J  =  1, NPUMP 
LINK  ( J ) =0 
00  CONTINUE 

DO  90  Jl  =  1, NVAR 
ERATE( J1 ) =0. 0 
0  CONTINUE 


ENTER  PUMP  INFORMATION 

1)  LINK  NUMBER  OF  PUMP 

2)  EFFICIENCY  INFORMATION 

a)  FLOW 

b)  EFF. 

WRITE(IOUT, 1007) 

DO  102  1=1, NPUMP 

READ(1,  1009)LINKd>,  XO,  YO,  XI,  Yl,  X2,  Y2 
CALL  SCURVE 

COEFFd,  1 )  =C0E1 
COEFFd,  2 )  =C0E2 
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COEFF(I,  3)=C0E3 


WRITE  ( IOUT,  1008 ) COEFF (1,1),  COEFF (1,2),  COEFF  (1,3) 

DO  111  11=1,0(1) 

IF(IPKIl)  .EQ.  LINK  ( I )  )  THEN 
LFLAG( I ) =11 
GO  TO  102 
ENDIF 

111  CONTINUE 
102  CONTINUE 


C  ENTER  INITIAL  DATA  ON  DECISION  VARIABLES 

C 

C  1)  INITIAL  GUESS 

C  2)  ELECTRICAL  RATE 

C  3)  LOWER  BOUND 

C  4)  UPPER  BOUND 


GUESS=0.  0 
RL0WBC=0. 0 
RUPPBC=0. 0 
WRITE(IOUT,  1010) 

DO  103  1=1, NVAR 
READ ( 1,  1012 )  GUESS,  ERATEd) 

XX ( I ) =GUESS 
BLVAR(I)=0. 

BUVAR ( I ) =NSTEP 

WRITEdOUT,  1011)1,  XXd),  BLVAR(I),  BUVAR(I) 
103  CONTINUE 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 


c 

c 

c 


ENTER  CONSTRAINT  INFORMATION 

1)  PRESSURE  CONSTRAINT  INFORMATION: 

NNODES  -  NUMBER  OF  NODES  TO  BE  CHECK  FOR  CONSTRAINT  VIOLATION 
PMIN  -  MINIMUM  ACCEPTABLE  PRESSURE  (PSI) 

PMAX  -  MAXIMUM  ALLOWABLE  PRESSUER  (PSI) 

NNNODE ( I )  -  NODE  NUMBER  OF  NODE  TO  BE  CHECKED 

2)  STORAGE  TANK  CONSTRAINT  INFORMATION: 

NCON  -  CONNECTING  PIPE  NUMBER  TO  TANK 
NNODE  -  NODE  NUMBER  OF  THE  TANK 
BOTEL  -  BOTTOM  ELEVATION  OF  THE  TANK 
TOPEL  -  TOP  ELEVATION  OF  THE  TANK 


READ(1, 1037)NN0DES,  PMIN, PMAX 
READ(1,  1038) (NNNODE(I), 1  =  1, NNODES) 

DO  220  1=1, NNODES 
DO  221  11=1,0(2) 

IFdNOdl)  .EQ.  NNNODEd)  )THEN 

Figure  B2.  (Cont’d). 
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NNNODE( I ) =11 
GO  TO  220 
ENDIF 

221  CONTINUE 
220  CONTINUE 


DO  104  1=1, NTANK 

READ ( 1,  1013)NN0DE(I),  NCON(I),  TOPEL(I),  BOTEL(I), DIAM, ENDEL(I), 
1  STOMIN  < I ) , PSTEP 

AREA ( I ) = ( 3. 1415927/4. ) «< DIAM*DIAM ) 

DO  105  J= 1, 0(2) 

IFCNNODE(I)  .  EQ.  INO(J))  THEN 
CTANK ( I ) =  D0 ( J ) 

NFLAG( I ) = J 
GO  TO  104 
ENDIF 

105  CONTINUE 
104  CONTINUE 


NN2=NN1 *NNODES 
NN3=NN1»NTANK 
L44  =  0 

DO  223  K=l, NN1 

DO  222  1=1, NNODES 
L44=L44+ 1 
BLCON ( L44 ) =PMIN 
BUCON (L44)=PMAX 

222  CONTINUE 

DO  224  11=1, NTANK 
L44=L44+1 

BLCON  ( L44 )  =STOMIN  (ID 
BUCON ( L44 ) =TOPEL < II ) *0. 55 
C  IF ( K  .EQ.  NN1 ) BLCON ( L44 ) =ENDEL( I 1 ) 

224  CONTINUE 

223  CONTINUE 


ENTER  INFORMATION  ON  WATER  USE  PATTERN.  ENTER  RATIO  OF  WATER  USE 
FOR  GLOBAL  DEMAND  MULTIPLICATION.  ENTER  ONE  VALUE  FOR  EACH  TIME 
C  STEP,  NN1 

READ(1,  1034) (RATIO(L),  L=l, NN1) 

C  ENTER  NODES  OF  WATER  USE  IN  AN  ARRAY  WHICH  WILL  NOT  CHANGE 


L45=0 

DO  99  11  =  1,0(2) 

IF( DO ( Ml )  .LE.  0.  .OR.  D0(M1)  . GT.  9E9)G0  TO  99 

Figure  B2.  (Cont’d). 
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non  non  no 


L45-L45+ 1 
CDEM ( L45) =D0( Ml ) 
99  CONTINUE 


C  DETERMINE  INTIAL  TANK  ELEVATIONS 


LTANK= 1 

WRITE ( IOUT, 1014) 
CALL  GC0MP(G, XX) 
LTANK=0 
ID0N=0 


0TPRNTr. FALSE. 
INPRNT=. TRUE. 
0TPRNT=. TRUE. 
NC0RE=20000 
NNVARS-NVAR 
NFUN=NN2+NN3+1 
MAXBAS=NFUN 
MAXHES=NVAR 
NNOBJ  =  NFUN 


READ  IN  TTITLE 
READ ( L, 1030 ) TTITLE 
WRITE ( IOUT, 1030 ) TTITLE 


ENABLE  ALL  DEFAULT  PARAMETERS 
C 

DO  106  K3=  1, 19 
DEFAUL ( K3 ) = 1. 0 
106  CONTINUE 

C  ENABLE  INPUT  FOR  TOLERANCE  ON  FRACTIONAL  CHANGE  IN  OBJECTIVE 
C  FUNCTION  (FPSTOP) 

C 

C  ENABLE  INPUT  FOR  CHANGE  IN  METHOD  OF  COMPUTING  PARTIAL 
C  DERIVATIVES  ( LDERIV :  O-FOWARD,  1-CENTRAL,  2-USER  SUPPLIED) 

DEFAUL( 1 ) =2.  0 
DEFAUL ( 2 ) =2.  0 
DEFAUL ( 3 ) =2.  0 
DEFAUL(6) -2.  0 
DEFAUL(7 ) -2.  0 
DEFAUL < 9) =2.  0 
DEFAULf 10 ) =2.  0 
DEFAUL( 15 ) =2.  0 
DEFAUL( 16) =2.  0 

Figure  B2.  (Cont  d). 
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DEFAUL ( 17  >  =2. 0 

READ< 1,  1031 )FPNEWT, FPINIT,  FPSTOP 

C  SET  TOLERENCE  PARAMETERS  TO  DEFAULT  VALUES  IF  READ  VALUES  OF  ZERO 

IF ( FPNEWT  .EQ.  0. ) DEFAUL ( 1 ) =1. 0 
IF (FPINIT  .EQ.  0. )DEFAUL<2)=1.0 
I F ( FPSTOP  .EQ.  0. )DEFAUL(3)=1.0 

READ ( 1, 1032 ) NNSTOP, IITLIM,  IIPR,  IIPN4,  IIQUAD,  LDERIV,  MMODCG 
C  SET  TOLERENCE  PARAMETERS  TO  DEFAULT  VALUES  IF  READ  VALUES  OF  ZERO 

I F ( NNSTOP  .EQ.  0. ) DEFAUL ( 6 ) =1. 0 
IF ( IITLIM  .EQ.  0. )DEFAUL(7)=1.0 
IF( IIPR  .EQ.  0. )DEFAUL<9)=1.0 
IF ( IIPN4  .EQ.  0. )DEFAUL(10)=1. 0 
IF ( IIQUAD  .EQ.  0.  )DEFAUL< 15 >  =  1. 0 
IF( LDERI V  .EQ.  0. ) DEFAUL ( 16 ) =1. 0 
IF ( MMODCG  .EQ.  0.  >DEFAUL( 17 ) =1. 0 

C  ECHO  VALUES 

WRITEdOUT,  1033)FPN£WT,  FPSTOP,  NNSTOP,  IITLIM,  IIPR,  IIPN4,  IIQUAD, 
1LDERIV, MMODCG 

C 

C 

C  CALL  OPTIMIZATION  ROUTINE 

C 

1000  CALL  GRGSUB ( INPRNT, OTPRNT, NCORE, NNVARS, NFUN, MAXBAS, MAXHES, NNOBJ, 
1TTITLE, BLVAR, BUVAR,  BLCON,  BUCON,  DEFAUL,  FPNEWT,  FPINIT,  FPSTOP, 
2FPSPIV, PPH1EP, NNSTOP, IITLIM, LLMSER, IIPR, IIPN4, IIPN5, IIPN6, IIPER, 
3IIDUMP, IIQUAD, LDERIV, MMODCG,  RAMCON,  RAMVAR,  XX,  FCNS,  INBIND, 

4RMULTS, NONBAS, REDGR,  NBIND, NNONB, INFORM, Z) 

C 

C 

C  DETERMINE  FINAL  TANK  ELEVATIONS 

LTANK=1 

WRITEdOUT,  1014) 

CALL  GCGMP(G, XX) 

C 

CALL  DP 

CALL  PRNOUT ( 'C'  ) 

C 

DO  998  L3=l, NVAR 

998  IF( XX ( L3)  .  LT.  0. 25 ) XX ( L3 ) =0.  0 
LTANK=2 

WRITEdOUT,  1014) 

CALL  GCOMP( G, XX) 

IBUG=0 

IF (LOOP  .LT.  LPFLG)GO  TO  1000 
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c  *#»»*»*#*#**#*  FORMAT  BLOCK  *****«****##«*» 

C 

1001  FORMAT (IX, 37(1H»), /, IX, 'ENTER  THE  NAME  OF  INPUT  FILE’,/) 

1002  FORMAT ( A14 ) 

1003  F0RMAT(1X,37!1H«>, /, IX, 'ENTER  THE  NAME  OF  OUTPUT  FILE',/) 

1004  F0RMAT(2X, 'THE  PROGRAM  CANNOT  ACCESS  THIS  FILE’) 

1005  FORMAT (4I5,2F10.2,2I5) 

1006  FORMATS///,  2X,  'TOTAL  NUMBER  OF  PUMPS  OPTIMIZED :',  14,  / 

12X,  'TIME  STEP  =  ',15,'  HOURS',/ 

22X, 'TOTAL  NUMBER  OF  DECISION  VARI ABLES : ' , 17, /, 

32X, 'NUMBER  OF  TANKS  CONTRIBUTING  TO  HEAD: ',17/) 

1007  FORMAT (///, 20X, 'PUMP  EFFICIENCY  INFORMATION', //, 8X, 'COEFFICIENT  #1 
1', 6X, 'COEFFICIENT  #2 ', 6X, ' COEFFICIENT  #3',//) 

1008  FORMAT ( 3F20. 8 ) 

1009  FORMAT ( 15,  6F10.  2 ) 

1010  FORMAT(///, 28X, 'PUMP  GROUP  INFORMATION' ,//, 2X, ' PUMP  GROUP', 10X, 

1 ' INITIAL  GUESS', 10X, 'LOWER  BOUND ', 10X, ' UPPER  BOUND’,/) 

1011  FORMAT (4X, 13, 14X, F10. 2,  12X, F10. 2,  11X, F10. 2) 

1012  FORMAT ( 2F10. 2 ) 

1013  F0RMAT(2I5, 6F10. 2) 

1014  FORMAT !  /// ) 

1030  FORMAT ( A90 ) 

1031  FORMAT ( 3F15. 14 ) 

1032  FORMAT ( 715 ) 

1033  FORMAT ( ///, 2X, 65( 1H* ) , /, 25X, 'OPTIONS  DISPLAY 2X,  'CONSTRAINT  TOL 
1ERANCE  =  ', F15. 14, /, 2X,  'OBJECTIVE  FUNCTION  TOLERANCE  =  \  F15. 14, /, 
22X, 'NUMBER  OF  OBJECTIVE  CHANGE  ITERATIONS  =  ’,15, /, 2X, 

3 'NUMBER  OF  NEWTON  ITERATIONS  =  ', 15, /, 2X, ' PRINT  OPTION  =  ',15,/, 
42X, 'GRADIENT  PRINT  OPTION  =  ’, 15, /, 2X, ' BASIC  VARIABLE  SEARCH  =  ', 
515, /, 2X, 'METHOD  OF  OBTAINING  PARTIAL  DERIVATIVE  =  ', 15, /, 2X, 

6 ' USE  OF  CONJUGATE  GRADIENT  =  ', 15, /, 2X, 65! 1H*  ),////// ) 

1034  FORMAT ( 8F10. 0 ) 

1035  FORMAT!//, IX, 37 ( 1H»), /, IX, 'ENTER  THE  NAME  OF  SYSTEM  DATA  FILE',/) 

1036  FORMAT! IX, 37(1H*>,  /,  IX, 'ENTER  THE  NAME  OF  INFO  DUMP  FILE',/) 

1037  FORMAT! 15, 2F10.2) 

1038  FORMAT! 1615) 

27  F0RMAT(32X, ' INITIAL  GCOMP  CALL') 

C 

C 

CLOSE! 1) 

CLOSE ( 8 ) 

CLOSE ( I IN ) 

CLOSE! IOUT) 

999  STOP 
END 


SUBROUTINE  SCURVE 

C 

c 

C  THIS  SUBROUTINE  FITS  A  QUADRATIC  CURVE  THROUGH 

C  USER  SUPPLIED  POINTS  USING  LAGRANGIAN  POLYNOMIALS 

C 

Figure  B2.  (Cont’d). 
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IMPLICIT  REAL*8(A-H, 0-Z) 
IMPLICIT  INTEGER*4( I-N) 

COMMON/ BLK3/X0,  YO,  XI, Yl, X2, Y2 
COMMON /BLK4 /COE 1, C0E2, C0E3 
COMMON/FILE/ I IN, IOUT 


C 

C 

XXO=(XO-Xl ) *<X0-X2) 

RO=YO/XXO 

XX1=(X1-X0)*(X1-X2) 

RI=Y1/XX1 

XX2= ( X2-X0 ) » ( X2-X1 ) 

R2=Y2/XX2 

C 

CGE1=R0»X1 *X2+R1 *X0*X2+R2*X0*X1 
C0E2= ( -RO» ( X1+X2) ) - <R1* (X0+X2) ) - ( R2* (XO+X1 ) ) 
C0E3=R0+R1 +R2 
C 

RETURN 

END 


Figure  B2.  (Cont’d). 
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APPENDIX  C: 


GCOMP  SUBROUTINE 


The  GCOMP  subroutine  computes  the  value  of  the  objective  function  and  problem  constraints.  This 
subroutine  is  required  by  GRG2.  The  array  containing  the  most  current  values  of  the  decision  variables, 
X,  is  passed  to  the  GCOMP  routine  for  computation  of  the  objective  function  and  constraints.  An  array 
containing  the  computed  values  of  the  objective  and  problem  constraints,  G,  is  passed  back  to  GRG2  along 
with  the  array  containing  the  current  decision  variables. 

Upon  entering  the  GCOMP  routine,  values  in  the  X  array  are  transferred  to  an  array  named  XXX. 
This  transfer  switches  entries  in  the  one-dimensional  X  array  with  NVAR  elements  to  a  two-dimensional 
array  also  with  NVAR  elements.  The  two-dimensional  array  is  an  NN1  x  NPUMP  array  where  NN1  is 
equal  to  the  number  of  time  intervals  in  the  24-hr  time  period  and  NPUMP  is  equal  to  the  number  of 
pumps  in  the  system.  NVAR  is  the  number  of  decision  variables  included  in  the  optimization. 

Essentially,  the  way  the  GCOMP  routine  works  is  that  for  each  time  interval,  the  pump  operating  the 
least  amount  of  time  is  selected  and  analyzed.  WADISO  is  called  and  the  system  is  balanced  under  the 
assumption  that  all  pumps  with  a  positive  value  in  the  XXX  array  are  operating.  If  WADISO  is  unable 
to  balance  the  network,  the  coefficient  array  used  by  WADISO  is  recomputed  and  the  network  is 
rebalanced. 

After  returning  to  GCOMP,  all  tanks  are  evaluated  to  determine  if  any  tanks  that  were  closed 
previously  should  be  opened.  A  tank  that  is  initially  closed  full  is  opened  when  the  head  at  the  beginning 
node  of  the  connecting  pipe  is  less  than  the  head  at  the  ending  node  (base  of  the  tank)  plus  the  water  level 
in  the  tank.  Conversely,  a  tank  that  is  initially  closed  empty  will  be  opened  if  the  head  at  the  beginning 
of  the  connecting  pipe  is  greater  than  the  head  at  the  end  of  the  connecting  pipe.  If  a  tank  is  opened, 
WADISO  is  called,  the  system  is  rebalanced,  and  the  tanks  arc  rcchecked  until  all  tanks  that  should  be 
open  arc  opened. 

After  all  tanks  are  checked  to  decide  if  any  should  be  opened,  a  check  is  made  to  determine  if  any 
should  be  closed.  This  is  accomplished  by  finding  the  amount  by  which  the  tank  drains  or  fills.  A 
separate  subroutine  called  UPDATE  computes  the  change  in  water  level  in  all  tanks  based  on  the  flow  rate 
in  tank  connecting  lines,  the  area  of  the  tank,  and  the  lime  of  analysis.  The  time  of  analysis  is  equal  to 
the  minimum  amount  of  time  a  pump  operates  during  the  lime  interval  evaluated.  For  example,  if  three 
pumps  operate  during  a  4-hr  lime  interval  with  pump  1  operating  for  1  hr,  pump  2  operating  for  2  hr,  and 
pump  3  operating  for  3  hr,  the  time  of  analysis  would  be  equal  to  1  hr. 

If  a  tank  docs  fill  or  drain  completely,  the  lime  to  closure  is  computed.  If  two  or  more  tanks  drain 
or  fill  completely,  the  time  to  closure  is  computed  for  each  lank  and  the  minimum  time  to  closure  is 
selected.  Tanks  arc  closed  by  closing  the  tank  connecting  line  and  setting  the  water  elevation  in  the  tank 
equal  to  the  bottom  or  top  elevation  of  the  tank.  The  water  levels  in  the  remaining  tanks  are  recomputed 
based  on  the  time  to  closure  if  only  one  tank  closes,  or  the  minimum  time  to  closure  if  multiple  tanks 
close.  The  time  of  analysis  is  then  set  equal  to  the  minimum  time  to  closure. 

The  cost  of  pumping  is  computed  in  subroutine  COSTC  called  from  GCOMP.  At  the  beginning  and 
end  of  a  time  interval,  the  kilowatt  consumption  of  each  pump  is  determined  based  on  the  pump  head  and 
pump  discharge  provided  by  WADISO,  and  the  pump  efficiency.  The  kilowatt  demand  for  each  pump 
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is  averaged  and  added  to  find  the  total  average  kilowatt  demand  for  the  pump  station  for  the  time  step 
used.  The  cost  of  pumping  is  found  by  multiplying  the  total  average  kilowatt  consumption  for  the  pump 
station  by  the  electrical  rate  in  effect  for  this  time  interval  and  the  time  of  analysis.  The  cost  of  pumping 
for  this  small  time  increment  is  added  to  the  cumulative  cost  of  pumping  thus  far  to  obtain  a  total  cost 
of  pumping  for  the  entire  24-hr  time  period. 

If  one  or  more  tanks  is  closed  during  the  time  step  analyzed,  then  the  time  of  analysis  is  set  equal 
to  the  minimum  time  to  closure.  If  this  is  the  case,  then  the  system  is  rebalanced  with  those  tanks  closed 
as  necessary  and  the  entire  process  of  tank  analysis  and  pumping  cost  analysis  is  repeated  until  the  end 
of  the  time  step  is  reached.  Once  this  point  is  reached  or  if  no  tanks  close  during  the  time  step,  then  the 
entry  in  the  XXX  array  just  analyzed  is  set  equal  to  zero. 

At  the  end  of  a  time  interval,  the  problem  constraints  are  computed.  Problem  constraints  are  the 
pressure  at  critical  nodes  and  the  water  level  in  each  tank.  Node  pressure  is  provided  by  WADISO, 
whereas  tank  water  level  is  obtained  from  the  UPDATE  subroutine.  At  this  point,  the  next  time  interval 
is  analyzed  and  the  entire  process  repeated  until  the  end  of  the  24-hr  time  period  is  reached. 

The  source  code  for  the  GCOMP  routine  is  shown  in  Figure  Cl.  The  FORTRAN  source  codes  for 
the  UPDATE  and  COSTC  subroutines  are  shown  in  Figures  C2  and  C3. 


SUBROUTINE  RCOMP(G, X) 


IMPLICIT  REAL*8( A-H, O-Z) 

IMPLICIT  INTEGER»4( I-N) 

INTEGER  PNL, 0 
REAL  INTLEV 

PARAMETER  (PNL=200, LNL=150, MNL=20, IA=2000> 

PARAMETER  (NDC=150, NTP=20, NPP=10) 

DIMENSION  G 1 151 ) , X ( 150 ) , XXX ( 24, 10), NCLOSE(NTP) 

DIMENSION  TIME(NTP), TLEV (NTP), TMPLEV ( NTP ) 

COMMON  /BLK1 /  NVAR,  NTANK,  NSTEP,  NNOBJ, Cl, NN1, NPUMP, CDEM( 100) 

COMMON  /BLK2/  ERATE ( NDC ) , AREA ( NTP ) ,  COEFF ( NPP, 3 ) , LINK ( NPP ) 

COMMON  /BLK2A/  CTANK(NTP),  NNODE(NTP),  NCON(NTP),  RATI0(24), 
lLFLAG(NPP), NFLAG(NTP) 

COMMON  /BLK5/  BLVAR < NDC > , BUVAR < NDC > , BLCON ( 2«NDC ) , BUCON < 2*NDC ) 
COMMON  /TANK/  TOPEL ( NTP ) ,  BOTEL ( NTP ) , ENDEL ( NTP ) , LTANK 
COMMON  /NODES/  0 ( 15 ) , HE ( PNL ) , D0( PNL ) , EL ( PNL ) , S ( PNL ) , G1 ( PNL ) 

COMMON  /PIPES/  A( IA), CP(PNL), DI(PNL), XL(PNL), HW(PNL) 

COMMON  /TOPOL/  I  BE (PNL), lEN(PNL), IPI(PNL), INO(PNL), IBI ( PNL ) , IEI < PN 
1L) 

COMMON  /PRINT/  IPM, IPP, IPE, NDFLG, IBUG 

COMMON  /JOB/  JOB 

COMMON  /FLAG/  R3M, R4M, IFOPP 

COMMON  /FILE/  IIN, IOUT 

COMMON  /BLK7/  NTNTNT 

COMMON  /BLK8/  NNODES, PMIN, PMAX, NNNODE ( 20 ) , NN2, NN3 

C  MltMOHIHMIMIIItMIHmOmlMMIIIMIIHM 

C 

C  LATEST  UPDATE  ON  THIS  FILE:  12  DEC  80 

C 

C  MUMIHimHnHIUMMIHHIHIIIlHHHIHHHI 


C  THIS  IS  A  SUBROUTINE  TO  CALL  A  NETWORK  SOLVING  ROUTINE  TO  EVALUATE 
C  THE  OPERATING  CONDITIONS  OF  A  DISTRIBUTION  NETWORK 


WRITE(»,»>'  INSIDE  GCOMP  ' 


C  XXX  IS  A  TWO  DIMENSIONAL  DUMMY  ARRAY  WITH  THE  FIRST  ELEMENT 
C  CORRESPONDING  TO  THE  CURRENT  TIME  INCREMENT  AND  THE  SECOND 

C  ELEMENT  CORRESPONDING  TO  THE  PUMP  OF  INTEREST 

C  SWITCH  X  ARRAY  TO  XXX  ARRY 


L3  =  0 

IF (NTNTNT  . EQ.  1)  GO  TO  204 
DO  200  L - 1 , NN1 
DO  20 t  L 1 = 1 , NPUMP 

Figure  Cl.  Source  code  for  the  GCOMP  routine. 
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L3=L3* 1 

IF(X(L3>  .LE.  BUVAR(L3) (THEN 
XXX(L, LI )  =X(L3) 

ELSE 

HOLD=X ( L3 ) 

XXX ( L, L1)=BUVAR(L3) 

IF ( L3*NPUMP  .LE.  NVAR )THEN 
XXX (L*l, LI ) =HOLD-XXX(L, LI ) +X (L3+NPUMP ) 

ELSE 

C  XXX(L, LI ) =BUVAR( L3) 

XXX(L, LI )=X(L3) 

ENDIF 

ENDIF 

IF(IBUG  . GT.  1 )WRITE(8, * ) 'TIME  STEP  PUMP  NO  DV’, L, LI, XXX(L, LI ) 
201  CONTINUE 
200  CONTINUE 
GO  TO  205 


C  TRANSFER  FOR  FINITE  DIFFERENCE  COMPONENT  OF  GCOMP 

204  L3=0 

DO  202  L= 1, NN1 
DO  203  Ll  =  l, NPUMP 
L3=L3* 1 

XXX(L,L1)=X<L3) 

203  CONTINUE 
202  CONTINUE 


C  REINITIALIZE  TANK  LEVELS  TO  SAME  LEVEL  AT  BEGINNING  OF  DAY 
C  OPEN  ALL  PIPES  CONNECTING  TANKS  AND  REINITIALIZE  TIME  AND 
C  TANK  LEVEL  ARRAYS 

205  DO  63  L2- l, NTANK 
J2=NFLAG( L2) 

D0( J2)=CTANK(L2) 

TEMPLEV= ( D0<  J2) / 1E10) - 100. 

TIME(L2) =0. 0 
TLEV  <  L2) -0. 0 
DO  64  13=1,0(1' 

I F ( IPI (L3)  .NE.  NCON ( L2 ) )G0  TO  64 
IF(TEMPLEV  .GE.  T0PEL(L2) )THEN 
CP(L3)=1E20 
NCL0SE(L2) =1 

ELSEIFI TEMPLEV  . LE.  B0TEL<L2) )THEN 
CP(L3)=1E20 
NCLOSE ( L2>  =2 
ELSE 

CP ( L3)  =  <  4.  72*XL( L3) ) / ( <  HW(L3) *  *  1. 852 ) *  <  DI ( L3) »*4. 87 ) ) 
NCLOSE (L2)=0 
ENDIF 

64  CONTINUE 
63  CONTINUE 


Figure  Cl.  (Cont’d). 
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o  o 


NTEN  IS  THE  MAXIMUM  NUMBER  OF  WADISO  CALLS  ALLOWED  PER  PUMP  TIME 
INCREMENT 

NTEN=NTANK»5 
I F ( NTANK  .LE.  2)NTEN=20 

C  INITIALIZE  COST  VARIABLES 

TCOST =0. 0 
COSTCO. 0 
L45  =  0 


C  PRIMARY  DO  LOOP 
DO  100  M  =  l, NN1 


C  INITIALIZE  DUMMY  AND  FLAG  VARIABLES 

LLFLG=0 
NTEST=0 
TEST  r0. 0 
TIME1=0. 0 
RTEST=0. 0 
NEND=0 
RK1=0. 0 
RK2=0. 0 
RKSUM=0. 0 

C  UPDATE  DEMAND  PATTERN  BY  MULTIPLYING  ALL  DEMANDS  BY  THE  GLOABAL 
C  DEMAND  FACTOR  (RATIO) 

L44  =  0 

IF ( RATIO*  M )  .EQ.  0.  ) RATIO ( M  )  =  1. 0 
DO  99  Min, 0(2) 

I F ( DO (Ml)  .LE.  0.  .OR.  D0(H1)  . GT.  9E9)G0  TO  99 
L44rL44*l 

D0(M1)=CDEM(L44)*RATI0(M) 

99  CONTINUE 


C  FIND  PUMP  WHICH  OPERATES  FOR  THE  SHORTEST  DURATION  THIS  TIME  STEP 

19?  CC^FLOAT ( NSTEP ) ♦ 100000. 

RMIN-CC 

DO  102  12-1, NPUMP 

IF ( XXX ( M, 12)  . LE.  0. )  THEN 
NTEST-NTEST* 1 
XL(LFL.AG(  12)  )  =  -  1 . 0 
GO  TO  102 

ELSEIF ( XXX ( M, 12 )  .GT.  0.  .AND.  XXX<H,I2>  . LT.  RMIN)  THEN 
IF(XXX(M, 12)  .EQ.  TEST)  THEN 

Figure  Cl.  (Cont’d). 
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XXX (M,  12 ) -0. 0 
NTEST^NTEST  *1 
XL(LFLAG(I2))=-1.0 
GO  TO  102 
ENDIF 

RMIN=XXX ( M,  12) 
MFLAG= 12 

XL(LFLAG( 12) )=0. 0 
ELSE 

XL(LFLAG(I2) )=0.0 
ENDIF 

102  CONTINUE 


C  ENABLE  FLAG  LLFLG  TO  INDICATE  NO  PUMPS  OPERATING  THIS  TIME  STEP 


IF ( NTEST  .  GE.  NPUMP)  LLFLG= 1 


C  RESET  COUNTER,  NCALL,  THE  NUMBER  OF  TIMES  WADISO  IS  CALLED 
C  PER  TIME  STEP 

NCALL=0 
TIMESUM-O. 0 

C  COMPUTE  TIME  INCREMENT  FOR  COST  COMPUTATION  AND  TANK  LEVEL  CHANGE 

98  IF (LLFLG  . EO.  1)  THEN 

TIME2-NSTEP 
ELSE 

TIME2=XXX ( M, MFLAG) 

ENDIF 

DELT=TIME2-TIME1 

CCC=RMIN-TIMESUM 

RMINN=CCC 


C  CHECK  TO  SEE  IF  NO  PUMPS  ARE  RUNNING  <LLFLG=1 ) 

IF( LLFLG  .EQ.  1)  THEN 
WRITE ( *, 2 ) M 

IF ( IBUG  .GT.  1 ) WRITE (8,  4 ) M 
IF( RTEST  .GE.  NSTEP)GO  TO  100 

C  ASSURE  CONNECTION  TO  FIXED  GRADE 

DO  95  KK: 1, NTANK 

IFUBUG  .GT.  1)WRITE(8,  *) 'TANK  NO.  NCLOSE ' ,  KK,  NCLOSE(  KK ) 
IF( NCLOSE ( KK )  .EG.  0  )G0  TO  94 

Figure  Cl.  (Cont’d). 


79 


IF( NCLOSE ( KK )  . EQ.  1JTHEN 
WRITE ( « , 1007) 

WRITE!*,  »> 'TANK  IS  FULL  AND  NO  PUMPS  ON:  OPENING  TANK  LINE' 
NCLOSE <KK>=0 
DO  54  14=1,0(1) 

IF ( IPI ( 14 )  .NE.  NCON(KK) )G0  TO  54 

CP(I4)  =  (4.72»XL(I4))/(HW(I4D*1.  852 *DI  ( 14) **4. 87) 

54  CONTINUE 
ENDIF 

95  CONTINUE 

ELSE 

WRITE (  »,  1)M,  LINK(MFLAG),  XXX (M,  rirLAG) 

I F ( I  BUG  .GT.  1 ) WRITE ( 8, 3)M, LINK<  MFLAG ) , XXX ( M, MFLAG ) 

ENDIF 

c  » »««»«.«»»•**»'»»***»»»• «»!»«' ««»«'*» »»*««*»»*»»*» «««ff »«*«««« 

94  ZZZ=XXX!M, MFLAG) 

DDD=0. 0 

WRITE!*,  *> 'BALANCE  BEFORE  TANK  LEVELS  ADJUSTED’ 

CALL  BALANCE ( NCALL, NTEN, M, ZZZ ) 

CALL  COSTCtDDD, RK1, RK2, RKSUM) 

93  IF(IBUG  .GT.  4)  CALL  PRNOUT('C') 

C  HHHMIHIIIMHIlUimHmMKimilHHHIHIMimimmi 


DMAXrO. 0 

C  OPEN  TANK  LINES 

DO  108  K=l, NTANK 
IF(NCLOSE(K)  .EG.  0)G0  TO  108 

DO  111  Nl=l, 0< 1 ) 

IFdPI(Nl)  .NE.  NCON( K )  )G0  TO  111 

K1=HFLAG(K) 

IBIDIBKNl ) 

IEI 1  =  IEI <  N1 ) 

T1=HE( IBI1 ) 

T2=HE( IEI1 ) 

IF ( NCLOSE ( K )  .EG.  1  .AND.  T1  . GT.  T2)G0  TO  108 
IF< NCLOSE! K)  . EQ.  2  .AND.  T1  . LT.  T2)G0  TO  108 

HDIFF=ABS(T1-T2) 

IF ( UDIFF  .GE.  DMAX)THEN 
DMAX=HDIFF 
K3  =  K 
K4  =  N 1 
T 1 1 =T1 
T21 -T2 


Figure  Cl.  (Cont’d). 
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IF(  I  BUG  .GT.  1  )  WRITE  (  8,  9502  )  K3,  K4,  IPI(K4),T11,T21,DMAX 
ELSE 
ENDIF 

111  CONTINUE 


108  CONTINUE 


IF(  I  BUG  .GT.  1  )  WRITE  (  8,  9503 )  K3,  NCLOSE  ( K3 ) , Til, T21,  DPI  AX 
IF ( DMAX  . EQ.  0. 0)G0  TO  305 

IF ( NCLOSE ( K3 )  . EQ.  1  .AND.  Til  . LT.  T21  .OR.  NCLOSE ( K3)  . EQ.  2 
1  .AND.  Til  .GT.  T21 )  THEN 

CP(K4)  =  (4. 72«XL!K4) ) / < HW< K4 > » *  1. 852»DI ( K4 ) * *4. 87) 

NCLOSE ( K  3  >  =  0 
I F ( I  BUG  . LE.  1 ) GO  TO  994 
WRITE! * , * ) 

WRITE! *,*) 'OPENING  TANK  NUMBER  ',K3 
WRITE<8,  »> 'OPENING  TANK  NUMBER  ',K3 
WRITE! •,  »  ) 

994  GO  TO  94 

ENDIF 


C  FIND  THE  'TANK  SLOPE'  FOR  EACH  STORAGE  TANK 

305  IF ( LLFLG  . EQ.  I )  THEN 
T2  =  NSTEP 
ELSE 

T2rXXX ( M, MFLAG ) 

ENDIF 

TINC-(T2-TEST)-TIME5UM 
IF! IBUG  . LE.  1 )  GO  TO  221 
WRITEffl.  * ) 

WRITE ( 8,  »  ) 

WRITE (8,  * ) ' T2  TEST  TIMESUM  TINC' 

WRITE (8,  *)T2, TEST, TIMESUM, TINC 
WRITE! 8,  »  ) 

WRITE! 0, • ) 

221  DO  106  K - 1 , NTANK 

CALL  UPDATEUINC,  K,  INTLEV,  WATLEV,  TTT  > 

TMPLEV ( K ) = I NTL EV 
TLEV ( K ) =WATLEV 
TIHE! K ) -TTT 

IF! TIME! K )  .LT.  0. 0 ) TIME ( K ) =0. 0 
106  CONTINUE 


C  FIND  THE  MINIMUM  TIME  A  TANK  VIOLATES  BOUNDARY 
DO  107  K= 1, NTANK 

IF! IBUG  .GT.  1 )WRITE(8,  * ) 'TANK  TLEV  TIME ', K, TLEV ( K ), TIME ( K ) 
IF ( NCLOSE ( K )  . NE.  0)G0  TO  107 

Figure  Cl.  (Cont’d). 
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IF ( TLEV ( K )  .GE.  BOTEL(K)  .AND.  TLEV(K)  .  LE.  TOPEL(K))GO  TO  107 
IF< TIME ( K )  .LT.  RMINN )  THEN 
RMINN-TIME( K ) 

NHOLD=K 

IF <1  BUG  . GT.  1 ) WRITE ( 8,  *> 'RMINN  NHOLD RMINN, NHOLD 
ENDIF 

107  CONTINUE 


C  TEST  TO  SEE  IF  ANY  TANK  LEVELS  WERE  VIOLATED  (RMINN  <  CCC) 
IF ( RMINN  .EQ.  CCC)GO  TO  105 


C  CLOSE  THE  TANK  LINE  WHICH  VIOLATES  A  BOUNDARY  FIRST 
DO  109  L= 1, 0< 1 ) 

IF(IPI(L)  .NE.  NCON ( NHOLD ) )  GO  TO  109 
CP ( L ) = 1 E20 
GO  TO  110 
109  CONTINUE 


C  ASSIGN  HEAD  VALUES  AT  TANK  NODES  AND  ENABLE  FLAG  INDICATING 

C  A  CLOSED  PIPE. 

C  NCLOSE = 1 ,  TANK  CLOSED  FULL 

C  NCLQSE=2,  TANK  CLOSED  EMPTY 

110  I =NHOLD 

J=NFLAG ( I ) 

TEMFLEV-TLEV ( I ) 

IFfTEMPLEV  .LE.  BOTEL ( I ) ) THEN 
WATLEV^BOTEL! I ) 

DO! J  >  =  <  BOTEL ( I)  +  100)»1E10 
HE( J ) = BOTEL! I ) *EL( J ) 

NCLOSE ( I ) -2 

WRITE ( »,  1009) NHOLD, TEMPLEV, BOTEL (I ),  <  < M- 1 ) *NSTEP > *RMINN 
C  WRITE ( 8, 1009)NH0LD, TEMPLEV, BOTEL(I), ( ( M - 1 ) *NSTEP ) ♦RMINN 
IF(IBUG  .GT.  1 )WRITE(8,  1009)NH0LD,  TEMPLEV, BOTELt I), 

1  < <M-1)*NSTEP>*TEST*RMINN 
ELSEIF ( TEMPLEV  . GE.  TOPEL ( I ) ) THEN 
WATLEV=TOPEL(I) 

D0(J)=(T0PEL(I)*100)»1E10 
HE ( J ) -TOPEL  < I ) »EL< J) 

NCLOSE (  I ) =1 

WRITE  (  «,  10 10)  NHOLD,  TEMPLEV,  TOPEL  ( I),  (  ( PI  - 1 )  »NSTEP  )  +RMINN 
C  WRITE! 8, 1010) NHOLD, TEMPLEV, TOPEL < I), ( ( M- 1 ) »NSTEP )♦ RMINN 
I F( I  BUG  .GT.  1)  WRITE (8,  1010) NHOLD, TEMPLEV, TOPEL! I ), 
l  < ( M- I ) *NSTEP) ♦ TEST* RMINN 
ENDIF 


IF! RMINN  .LE.  CCC ) DELT-RMINN 

Figure  Cl.  (Conl’d). 
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IF*  TINC  .LT.  RNINN)DELT-TINC 


DO  103  K=l, NTANK 
J2=NFLAG( K ) 

IF ( NCLOSE ( K )  . EO.  2)  THEN 
DO ( J2  >  =  ( BOTEL( K)*100)»1E10 
ELSEIF ( NCLOSE ( K )  . EO.  1)  THEN 
DO*  J2 )  =  ( TOPEL (K)*100)*1E10 
ELSE 

D0(J2)Z(TMPLEV(K)*100)»1E10 

ENDIF 

103  CALL  UPDATE*  RMINN,  K,  INTLEV, WATLEV, TTT) 
CALL  BALANCE ( NCALL, NTEN, M, ZZZ) 

CALL  COSTC*RMINN, RK1, RK2, RKSUM) 


C  FIND  COST  OF  OPERATING  PUMPS  FOR  TIME=DELT 

105  IF( RMINN  . LE.  CCC ) DELT-RM INN 
IF ( TINC  .LT.  RMINN)DELT=TINC 

C  TIME  ADJUSTMENT 

TIMESUM=TIMESUM*DELT 
IF(IBUG  . GT.  1)  THEN 
WRITE*  8,  * ) 

WRITE ( 8,  » ) ' TIME2  TIMESUM ' ,  TIME2,  TIMESUM 
WRITEO,  »  ) 

ENDIF 

AABS=ABS*TIME2 -TIMESUM) 

IF ( A ABS  .LE.  0.0001)  THEN 
TIME1 =TIME2 
ELSE 

TIME1 -RMINN 
ENDIF 


IF*  RMINN  .  LT.  CCC)GO  TO  98 


WRITE* *,«) 'BALANCE  AFTER  TANK  LEVELS  ADJUSTED' 

CALL  BALANCEtHCALL, NTEN, M, ZZZ) 

CALL  COSTCtDELT, RK1, RK2, RKSUM) 

IF(XXX(M, MFLAG)  . EO.  NSTEP  .OR.  LLFLG  .EG.  1)NEND=1 


IF*  NEND  . NE.  1 >G0  TO  332 
C  EVALUATE  PRESSURE  CONSTRAINTS 
DO  331  J=l, NNODES 


Figure  Cl.  (Cont’d). 
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L45  =  L45d 
J 1 =NNNODE ( J ) 

G(L45)=(HE(Jl)-ABS(EL(Jl)))/2. 31 
IF( G ( L45 )  .LT.  PMIN ) WRITE< »,  1038 ) INCH  J1 ) , G< L45 ) 
IF ( G ( L45 )  .GT.  PMAX)WRITE<  «,  1039)  IN0(  J1 ),  G(L45) 
331  CONTINUE 


C  PRINT  TANK  LEVEL  AT  EVERY  TIME  STEP 

332  IF( LTANK  . GE.  1  )THEN 
IF ( LLFLG  .EQ.  1 ) THEN 
WRITEdOUT,  2)M 
ELSE 

WRITEdOUT,  1)M,  LINK (MFLAG),  XXX (M,  MFLAG) 
ENDIF 

WRITEdOUT,  1035) 

ENDIF 

DO  113  1=1, NTANK 

IF ( NEHD  .EQ.  1 ) THEN 
L45  =  L45d 
G ( L45 ) =TLEV ( I ) 

ENDIF 

IF ( LTANK  .GE.  1 ) WRITE (IOUT,  1037 ) I, TLEV ( I ) 
113  CONTINUE 


C  COMPUTE  COST  OF  OPERATION  THIS  TIME  INCREMENT 
I F ( NEND  . HE.  1 )G0  TO  333 

TCOST  =  TCOST ♦ ( RKSUM* ( ERATE< ( ( M- 1 ) »NPUMP ) d ) > / 100.  ) 
IFdBUG  .GT.  1  )WRITE(8,  1002)TC0ST 


333  IF< LTANK  .  GE.  1  IWRITEdOUT,  2000) 

IF (LLFLG  . EQ.  1 )G0  TO  100 

RTEST  =  XXX <  M, MFLAG) 

TIME1=XXX(M, MFLAG) 

XXX ( M, MFLAG ) =0. 

TEST  =RMIN 
NTEST =0 

GO  TO  199 

100  CONTINUE 


G ( NNOBJ ) =TCOST 

WRITE( », 1008 )TCOST, PENALTY, G( NNOBJ ) 

IF( LTANK  .EQ.  2)WRITE( IOUT, 1008)TC0ST, PENALTY, G(NNOBJ) 

Figure  Cl.  (Corn’d). 
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IF  1 1  BUG  .GT.  1 )  WRITE! 8,  1008 )TCOST, PENALTY, G( NNOBJ ) 


IBUG= 1 

IF< I  BUG  .GT.  3 )STOP 


C  uiMimiuHi  FORMAT  BLOCK  »***«»*#»*»»*» 


1  FORMAT! /,2X, 'CYCLE  NUMBER ’,  13,  10X,  ' PUMP  NUMBER’, 15, 5X, 'TIME  THIS  P 
1UHP\  F10.  5) 

2  FORMAT!/, 2X, 'CYCLE  NUMBER', 13, 10X, 'NO  PUMPS  OPERATING') 

3  FORMAT!///, 2X, 'CYCLE  NUMBER ', 13, 10X, ' PUMP  NUMBER ', 15, 5X, ’ TIME  THIS 
1PUMP’, F10. 5, //) 

4  FORMAT!///, 2X, 'CYCLE  NUMBER ', 13, 10X, ' NO  PUMPS  OPERATING',//) 

1001  FORMAT!  2X,  '  ***«»<"""•  ERROR  »**»**  DUMMAY  ARRAY  NOT  SWITCHED') 

1002  FORMAT! //,2X,  'TOTAL  COST  OF  PUMPING  THIS  TIME  INCREMENT :', F10.  2 ) 

1004  FORMAT ( 2X, 110, 1SX, 110, 21X,  F10. 2) 

1005  FORMAT! //, 33X, 'FLOW  INFORMATION ',//, 4X, ' BEGINNING  NODE ' , 17X, ' ENDIN 
1G  NODE', 19X, 'FLOWRATE', /) 

1006  F0RMAT(2X, 'INITIAL  WATER  LEVEL:', 

1F10. 4, /, 2X, 'TANK  FLOW F10. 4, /, 2X, ’ FINAL  WATER  LEVEL : ' , F10. 4, /, 
22X, 'TOTAL  HEAD  AT  TANK ’ , F10. 4, // > 

1007  FORMAT!/,  15X,  '»*»*  WARNING  »»»»',/, 4X,  'NO  CONNECTION  TO  SUPPLY  P 
10INT,  MUST  TURN  OH  PUMP',/) 

1008  FORMAT <2X, 'COST  OF  PUMPING  THIS  24  HR.  PERIOD', F10. 2, /, 

12X, 'PENALTY  FUNCTION  THIS  24  HR.  PERIOD' , F10. 2, /, 

22X, 'OBJECTIVE  FUNCTION  THIS  TIME  PERIOD ’, FID. 2, III > 

1009  FORMAT! /, 20X, '**»«*»»**»»*»**  WARNING  *«»»»»«»*»»*»»»',/, 2X, 

1 'TANK  NUMBER', 14, 2X, 'HAS  ELEVATION  OF ' , F7. 3, 2X, ' WHICH  VIOLATES  LOW 
2ER  BOUND  CONSTRA I  NT ' , F7. 3, 2X,  'OCCURING  AT  TIME  =',F7. 4, /> 

1010  FORMAT!/, 20X, '**«»**»***»*»»»  WARNING  ***»*»»«*»»»»»«',/, 2X, 

1 'TANK  NUMBER', 14, 2X. 'HAS  ELEVATION  OF ' , F7. 3, 2X, ’ WHICH  VIOLATES  UPP 
2ER  BOUND  CONSTRAINT ', F7. 3, 2X, 'OCCURING  AT  TIME  =',F7. 4, /> 

1035  FORMAT!/, 10X, 10(1H*),3X, 'FINAL  TANK  WATER  LEVELS’ , 3X, 10! 1H* ), //, 
114X, 'TANK  NUMBER', 20X, 'ELEVATION  (FT)'/) 

1037  FORMAT! 15X, 15, 24X, F10. 3) 

1038  FORMAT! /, 20X,  '  ***«*****«»**»»  WARNING  »**«*»»»»«»»*»»',/, 2X, 

1 ’ NODE  NUMBER’, 15, 3X, 'HAS  VIOLATED  LOWER  PRESSURE  CONSTRAINT’,/, 

22X,  'PRESSURE  AT  THIS  NODE  =  \F10.  2, /> 

1039  FORMAT!/, 20X,  WARNING  »«»»«»»*»»*«*»»',/, 2X, 

1 ' NODE  NUMBER', 15, 3X, 'HAS  VIOLATED  UPPER  PRESSURE  CONSTRAINT',/, 

22X, 'PRESSURE  AT  THIS  NODE  =  ',F10. 2, /) 

2000  FORMAT!///) 

2001  FORMAT!////) 

2002  FORMAT!//, 2X, 'SIMULATION  AT  TIME  = ' , F10. 4, 3X, 'TANK  LINE  CLOSED') 

2003  F0RMAT(2X, 'TANK  LINE  IS  NOW  OPENED') 

3005  FORMAT!/, 2X, 'HEAD  AT  BEGINNING  NODE ’ , 15, 3X, F10. 3, /, 2X, ’ HEAD  AT  END 
1ING  NODE',  3X,  15,  3X,  F10.  3, /) 

9500  FORMAT!//, 2X, 'TESTING  NCLOSE  TANK  NUMBER', 15, 4X, 'NCLOSE', 15, //) 

9501  FORMAT!//, 2X, 'TANK  NO.  NCLOSE  B  NODE  E  NODE  T1  T2', 

1415, 2F10. 4) 

9502  FORMAT! 2X,  'K1  K2  LINK  NO.  Til  T21  DMAX' , 315, 3F10. 2) 

Figured.  (Corn’d). 
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9503  FORMAT* //, 2X, 'TANK  NO.  NCLOSE  T1  T2  DMAX ' , 215, 3F10. 4 ) 

RETURN 

END 


SUBROUTINE  CHECK 
WRITE*  »,  »  ) 

WRITE*  », »  ) 

WRITE* »,») 'TERMINATING  PROGRAM  DUE  TO  INSTABILITIES  IN  TANK  OPERAT 
1  ION ' 

WRITE*  »,  «  ) 

WRITE*  »,  «  ) 

WRITE  *  8,  «  ) 

WRITE(8,  «  ) 

WRITE(8,  •> 'TERMINATING  PROGRAM  DUE  TO  INSTABILITIES  IN  TANK  OPERAT 
1  ION ' 

WRITE*  8,  »  ) 

WRITE ( 8,  •  ) 

STOP 

RETURN 

END 


SUBROUTINE  BALANCE ( NCALL,  NTEN, M, ZZZ> 

IMPLICIT  REAL*8( A-H, 0-Z) 

IMPLICIT  INTEGER»4< I-N) 

INTEGER  PNL,  0 

PARAMETER  <PNL=200,  LNL=150,  MNL=20,  IA  =  2000) 

PARAMETER  (NDC=150,  NTP  =  20, NPP=10) 

COMMON  /BLK2/  ERATE ( NDC ) ,  AREA ( NTP ) , COEFF ( NPP, 3 ) , LINK ( NPP ) 

COMMON  /NODES/  0*15 ) ,  HE* PNL ) ,  DO* PNL ) ,  EL ( PNL) ,  S ( PNL ) ,  G1  ( PNL ) 

COMMON  /PIPES/  A* IA),  CP (PNL),  DI(PNL),  XL (PNL),  HW(PNL) 

COMMON  /TOPOL/  I  BE  < PNL),  IEN(PNL),  IPI(PNL),  INO(PNL),  IBI(PNL),  IEI(PN 
1L) 

COMMON  /PRINT/  IPM, IPP, IPE, NDFLG, IBUG 
COMMON  / ACCU/  PRAC, FLAC, HWMA, ICL 
COMMON  /JOB/  JOB 
COMMON  /FLAG/  R3M, R4M, IFOPP 
COMMON  /FILE/  I  IN,  IOUT 

ISUMrO 


C  CALL  SIMULATION  ROUTINE  TO  BALANCE  THE  NETWORK 


Figured.  (Cont’d). 
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u  o  cr»  u  u  o  u 


WRITE<»,  M'CALLING  WADISO' 

CALL  SIMBALt 'O' ) 

ISUM= ISUM* 1 

C  CHECK  TO  SEE  IF  EXCEED  MAXIMUM  NUMBER  OF  CALLS 
NCALL=NCALL*1 

IF( NCALL  .GE.  NTEN )  CALL  CHECK 


CHECK  TO  ASSURE  SYSTEM  IS  BALANCED,  IF  NOT  RECONSTRUCT  SPRASE 
MATRIX  AND  REBALANCE. 

IF ( 0 ( 9 )  . EQ.  3)  RETURN 

IFOPP  =  100,  WADISO  FAILURE  DUE  TO  INFINITE  LOOP 

IFOPP  =  997,  WADISO  FAILURE  DUE  TO  DIVISION  BY  ZERO 

IFOPP  =  998,  WADISO  FAILURE  DUE  TO  INABILITY  TO  DELIVER  HEAD 

IFOPP  =  999,  WADISO  FAILURE  DUE  TO  INABILITY  TO  DELIVER  DISCHARGE 

IF ( IFOPP  .EO.  996)  GO  TO  88 
IF ( IFOPP  .EO.  997)  GO  TO  84 
IF ( IFOPP  .EO.  998)  GO  TO  84 
IF ( IFOPP  .EQ.  999)  THEN 
RETURN 
ENDIF 

IF( IFOPP  .EQ.  100)  GO  TO  84 

88  WRITE* », 1011 ) 

IF( IFOPP  .NE.  998)0(4 ) =0 

CALL  SIMBAL( '0' ) 

ISUM=ISUM* 1 
IF ( 0 ( 9 )  .EQ.  3)  THEN 
WRITE( »,») 'SUCCESSFUL  REBALANCE' 

RETURN 

ENDIF 

C  ONE  LAST  TRY  TO  REBALANCE 

IF  ( ISUM  .GT.  2 )G0  TO  90 
84  0( 4 ) =0 

IFUBUG  .GT.  1)WRITE(8,  «) 'ISUM  IFOPP',  ISUM,  IFOPP 
DO  91  K9K= 1,0(1) 

IF(CP(K9K)  .GT.  0)  THEN 
A ( K9K ) = 1 , /CP ( K9K ) 

ELSE 

Y 1 =CP ( K9K ) 

Y2rDI ( K9K ) 

Y3=HW ( K9K ) 

A ( K9K ) = 1 . 85 /SORT ( ( Y2* Y2 ) -2. »Y1»Y3) 

ENDIF 

91  CONTINUE 


Figure  Cl.  (Cont’d). 
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GO  TO  80 

90  WRITE(8, 1 )  M, LINK(MFLAG), ZZZ 

ITMPr ICL 
ICLr 100 

CALL  SIMBAL! *0’ ) 

ICL^ITMP 

I F <  0 ( 9 )  .EQ.  3 ) THEN 

WRITE! », ■) 'SUCCESSFUL  REBALANCE' 

RETURN 

ENDIF 

WRITE ( 8,  ») 

IF ( I  BUG  .GT.  1)CALL  PRNOUT('C') 

WRITE! 8, 5IR3N/2. 308, R4M»448. 8 
WRITE!  »,  5JR3H/2.308,  R4HM4Q.8 

C  . . .  FORMAT  BLOCK  »»»**• *»»»»** **««»« 

1  F0RMAT(2X, 'CYCLE  NUMBER ',  13,  10X,  ' PUMP  NUMBER ', 15, 5X, ' TIME  THIS  PUM 
IP’, F10. 5) 

5  FORMAT! //, 2X, 'MAXIMUM  PRESSURE  IMBALANCE  =  ', F6, 2, 4X, ' MAXIMUM  FLOW 
1  IMBALANCE  =  ',F6. 2, //) 

10U  FORMAT!/, 2X, 'RECONSTRUCTING  SPARSE  MATRIX  AND  REBALANCING’,/) 


RETURN 

END 


Figured.  (Cont’d). 
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SUBROUTINE  UPDATE ( DELT, I,  INTLEV,  WATLEV,  TTT ) 

IMPLICIT  REAL»8( A-H,  O-Z) 

IMPLICIT  INTEGER*4 ( I -N ) 

INTEGER  PNL,  0 
REAL  INTLEV 

PARAMETER  vPNL=200,  LNL=150,  MNL=20, IA  =  2000) 

PARAMETER  ( NDC  = 150, NTP=20, NPP=10) 

COMMON  /BLK1/  NVAR,  NTANK,  NSTEP,  NNOBJ,  Cl, NN1, NPUMP, CDEM ( 100 > 

COMMON  /BLK2/  ERATE<  NDC5 ,  AREA <  NTP ),  COEFF < NPP,  3) ,  LINK( NPP ) 

COMMON  /BLK2A/  CTANK(NTP), NNODE(NTP), NCON(NTP), RATIO (24), 
lLFLAG(NPP),  NFLAG(NTP) 

COMMON  /BLK5/  BLVAR(NDC), BUVAR(NDC), BLCON ( 2*NTP ) ,  BUC0N(2»NTP) 

COMMON  /TANK/  TOPEL ( NTP ) , B0TEL<  NTP ) , ENDEL, LTANK 

COMMON  /NODES/  0< 15 ) ,  HE ( PNL ) ,  DO < PNL ) , EL ( PNL ) , S< PNL ) , G1 ( PNL ) 

COMMON  /PIPES/  A( IA>, CP(PNL), DI(PNL), XL(PNL), HW(PNL) 

COMMON  /TOPOL/  IBE(PNL),  IEN(PNL>,  IPI(PNL), INO(PNL), I BI ( PNL ) , IEI ( PN 
1L) 

COMMON  /PRINT/  IPM, IPP, IPE, NDFLG, IBUG 
COMMON  /JOB/  JOB 
COMMON  /FLAG/  IFOPP 
COMMON  /FILE/  IIN, IOUT 


J=NFLAG( I ) 

C  COMPUTE  TANK  LEVEL  FROM  WADISO  ARRAY  DO 
WATLEV  =(D0(J)/1E10)-100. 

C  FIND  THE  FLOW  IN  THE  TANK  CONNECTING  LINE 
DO  106  K-l, 0( 1 ) 

IF(IPI(K)  .NE.  NCON ( I ) )G0  TO  106 

C  FIND  HYDRAULIC  GRADE  AT  BEGINNING  AND  ENDING  NODES 

IBI 1  =  1 BI <  K ) 

IEI 1  =  IEI <  K ) 

T1 =HE  C IBI 1 ) 

T2=HE ( IEI 1 ) 

DELH-ABS(T1-T2) 


C  COMPUTE  FLOWRATE 

TFLOW= ( DELH/CP< K ) ) » »0.  54 
FILL= ( 3600. * TFLOW* DELT) /AREA ( I ) 

INTLEV=WATLEV 

Figure  C2.  Source  code  for  the  UPDATE  routine. 
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c 

c 

c 

106 

C 

998 


C 

1001 

1002 

1003 

1004 

1005 


DETERMINE  DIRECTION  OF  FLOW,  NOTE  TANK  IS  ALWAYS  ENDING  NODE 

IF(T1  .GT.  T2)  THEN 

TANK  IS  FILLING 
C33  1. 

ELSE 

TANK  IS  DRAINING 
C3--1. 

ENDIF 
GO  TO  998 
CONTINUE 

UPDATE  TANK  LEVEL:  WATLEV 

C4=B0TEL( I) 

C5=T0PEL ( I ) 

WATLEV3 INTLEV* <C3*FILL> 

D0(J)=(WATLEV*100)»1E10 
HEADLEV=WATLEV*EL< J> 

IF ( DELT  .GT.  0)  THEN 
SLOPE3 (C3* FILL) /DELT 
IF <  SLOPE  .GE.  0)THEN 
TTT3 ( C5- INTLEV ) /SLOPE 
ELSE 

TTT= ( INTLEV-C4) /ABS( SLOPE) 

ENDIF 
ELSE 
TTT=0. 0 
ENDIF 


IF(IBUG  .GT.  1.  AND.  I  . LE.  1 )WRITE<8, 1004 >DELT 
IF( I  BUG  .GT.  1)WRITE<8, 1005)1, TFL0W*448.  84,  INTLEV,  <C3»FILL), 
1  WATLEV, HEADLEV 


FORMAT  BLOCK  »«*»*«#»«»**** 

FORMAT ( 216, F13. 4, 16, 6F10. 4) 

FORMAT< // /, '  TANK  NO.  B  NODE  HEAD  E  NODE  HEAD 

1',  /) 

FORMAT (16,  110,  F14. 4, 6X, 16, FI 4. 4) 

FORMATS,  5X,  'TIME  STEP  THIS  PERIOD ',  F10.  4, 

1//,  '  TANK  NO.  FLOWRATE  INTLEV  DEL  H  WATLEV 

2  HEAD  AT  TANK', /) 

FORMAT ( 16,  5F13.  4) 

RETURN 

END 


Figure  C2.  (Cont’d). 


90 


SUBROUTINE  COSTC(DELT, RK1, RK2, RKSUM) 

IMPLICIT  REAL*8<  A-H, O-Z) 

IMPLICIT  INTEGER»4 ( I  - N  > 

INTEGER  PNL, 0 

PARAMETER  <PNL=200,  LNL=150,  MNL=20, IA=2000) 

PARAMETER  (NDC=150,  NTP=20, NPP=10> 

COMMON  /BLK1/  NVAR,  NTANK,  NSTEP,  NNOBJ,  Cl, NN1, NPUMP, CDEM( 100) 

COMMON  /BLK2/  ERATE < NOC > , AREA ( NTP > ,  COEFF < NPP, 3 ) , LINK ( NPP ) 

COMMON  /BLK2A/  CTANK ( NTP ) , NNODE ( NTP ) , NCON ( NTP > ,  RATIO< 24 ) , 
lLFLAG(NPP),  NFLAG(NTP) 

COMMON  /BLK5/  BLVAR < NDC ) , BUVAR ( NDC ) ,  BLCON ( 2*NTP ) ,  BUCON ( 2»NTP ) 

COMMON  /TANK/  TOPEL ( NTP > , BOTEL < NTP > ,  ENDEL,  LTANK 

COMMON  /NODES/  0(15 ) ,  HE ( PNL ) ,  D0< PNL ) ,  EL< PNL ) , S ( PNL) ,  G1 (PNL ) 

COMMON  /PIPES/  A( I A ) ,  CP(  PNL) ,  DI  ( PNL) ,  XL(PNL),  HVI  ( PNL) 

COMMON  /TOPOL/  IBE(PNL),  IEN(PNL),  IPI(PNL),  INO(PNL), IBI(PNL), IEI(PN 
1L) 

COMMON  /PRINT/  IPM,  IPP, IPE, NDFLG, I  BUG 
COMMON  /JOB/  JOB 
COMMON  /FLAG/  IFOPP 
COMMON  /FILE/  IIN,IOUT 


I F ( RK2  . LE.  0. 0)G0  TO  9998 
RK1 -RK2 


C  FIND  FLOW  AND  HEAD  ASSOCIATED  WITH  EACH  PUMP  THAT  IS  OPERATING 

9998  IF(IBUG  .  GT.  1 )WRITE<8, 1012) 

RKW=0. 0 

DO  103  1=1, NPUMP 
J1=LFLAG( I ) 

IF( A ( J1 )  , EQ.  0. )  THEN 
QADD=0. 

ELSE 

71=CP( J1 ) 

Z2=DI(J1) 

QPCFS=(-Z2-1.85/A< Jl))/2/Zl 
QFGPM  =  QPCFS»448.  8 

EFFI =COFFF ( I, 1 ) ♦ <COEFF( I,  2) »QPGPN ) ♦ ( COEFF (1,3) «QPGPM»QPGPM ) 
OADD  =  QPCFS/ ( EFFI/100.  ) 

ENDIF 


C  FIND  PUMP  HEAD 

HI =HE( I B I ( J 1 ) ) 

H2  =  HE< IEI <  J 1 ) ) 

HFAD=ABS(H1 -H2) 

Figure  C3.  Source  code  for  COSTC  routine. 
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C  COMPUTE  KILOWATT  DEMAND  FOR  EACH  PUMP  AND  TOTAL  PUMP  STATION 

C  KILOWATT  DEMAND 

ZKW=Cl»HEAD»OADD 

RKW=RKW*ZKW 

IF ( IBUG  .GT.  1  .AND.  ZKW  . GT.  0.  > WRITE ( 8, 101 1 > LINK < I > , HEAD, 

1  QPGPM, EFFI, ZKW 

103  CONTINUE 

RK2-RKW 

IF ( RK 1  . LE.  0.0  .AND.  DELT  .  LE.  0.  0 ) RKAREA-O.  0 
IF ( RK 1  . LE.  0. 0)G0  TO  9999 
RKAREA= ( DELT* (RK1-RK2) 1/2. 

RKSUM=RKSUM*RKAREA 

9999  IF(IBUG  . GT.  I )WRITE*8, 1002) 

I F ( I  BUG  .GT.  1 ) WRITE (8,  1003) HEAD,  RK1, RK2, DELT, RKAREA 
IF(IBUG  .GT.  1 >WRITF<8, 1004)RKSUH 

C  . . .  FORMAT  BLOCK  *»*»»»*****»»» 

1002  FORMAT ( //, 33X,  'PUMP  INFORMATION //9X,  ' HEAD ’, 9X,  ' KW  #1  USED',4X, 
l’KW  #2  USED’, 7X,  'DELTA  T ' , 7X,  ' AVG.  KW  USED',/) 

1003  FORMAT ( 5F15. 4 ) 

1004  FORMAT*/, 5X, 'ACCUMULATED  KILOWATT  USAGE : ' , F10. 4 ) 

1011  FORMAT* 2X,  15,  11X,  F10. 2, 10X,  F10. 2, 5X, F10. 2, 4X, F10. 2) 

1012  FORMAT*//, 2X, 'PUMP  NUMBER 8X, ' PUMP  HEAD ', 7X, ' PUMP  DISCHARGE’, 
15X, 'EFFICIENCY', SX, 'KW', //) 


RETURN 

END 


Figure  C3.  (Cont’d). 
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APPENDIX  D 


FORT  HOOD  HYDRAULIC  NETWORK  DATA 

Table  D1 

Pipe  Data  for  Fort  Hood  Main  Post 


Pipe 

Number 

Connecting  Nodes 

Diameter 

(In.) 

Length 

(ft) 

Hazen-VVIlliams 

C-Factor 

1 

1 

2 

30 

200 

100 

2 

2 

3 

30 

1 

100 

3 

2 

5 

30 

1 

100 

4 

2 

7 

30 

1 

100 

5 

2 

9 

30 

1 

100 

6 

2 

11 

30 

1 

100 

7 

3 

4 

Pump 

- 

- 

8 

5 

6 

Pump 

- 

- 

9 

7 

8 

Pump 

- 

- 

10 

9 

10 

Pump 

- 

- 

11 

11 

12 

Pump 

- 

- 

12 

4 

13 

30 

1 

100 

13 

6 

13 

30 

1 

100 

14 

8 

13 

30 

1 

100 

15 

10 

13 

30 

1 

100 

16 

12 

13 

30 

1 

100 

20 

13 

15 

30 

200 

100 

21 

15 

16 

24 

120 

100 

22 

16 

18 

24 

200 

100 

2T 

15 

17 

30 

160 

100 

24 

17 

18 

30 

195 

100 

25 

16 

19 

12 

2800 

35 

26 

19 

20 

12 

1020 

35 

27 

20 

25 

12 

100 

100 

28 

20 

21 

8 

1350 

35 

29 

19 

70 

8 

880 

35 

30 

21 

22 

8 

3220 

35 

31 

22 

23 

12 

1500 

35 

32 

23 

90 

12 

5060 

35 

33 

18 

24 

24 

2520 

84 

34 

24 

65 

24 

1800 

84 

35 

18 

26 

16 

210 

69 

36 

26 

27 

8 

3860 

100 

37 

26 

28 

16 

2080 

69 

38 

28 

72 

16 

900 

69 

39 

29 

31 

16 

3080 

69 

41 

31 

32 

16 

50 

100 

42 

33 

36 

8 

9733 

100 

43 

31 

34 

20 

9420 

96 

44 

34 

35 

20 

3820 

96 

46 

37 

75 

24 

3100 

84 

47 

37 

40 

18 

1020 

102 

48 

36 

74 

12 

2130 

102 

49 

38 

39 

12 

135 

100 

93 

Table  D1  (Cont’d) 


Pipe 

Number 

Connecting  Nodes 

D!  imeter 
(in.) 

Length 

(ft) 

Hazen-Wlllianis 

C-Factor 

50 

38 

40 

16 

520 

102 

51 

40 

42 

12 

4360 

102 

52 

37 

76 

18 

2275 

112 

53 

41 

42 

10 

1600 

112 

54 

41 

43 

16 

1840 

112 

55 

43 

46 

12 

2280 

130 

56 

17 

44 

20 

9372 

130 

57 

44 

45 

20 

60 

100 

58 

44 

46 

12 

2070 

130 

54 

46 

79 

10 

1960 

130 

60 

47 

48 

10 

1160 

130 

61 

43 

78 

16 

2090 

130 

62 

48 

49 

10 

3130 

130 

63 

44 

50 

18 

80 

100 

64 

47 

51 

8 

3640 

103 

65 

49 

52 

18 

540 

130 

66 

52 

53 

12 

1620 

130 

67 

53 

54 

8 

1020 

130 

68 

41 

53 

8 

2440 

130 

64 

42 

54 

8 

2480 

130 

70 

53 

55 

12 

5020 

130 

71 

54 

55 

8 

6230 

130 

72 

51 

56 

8 

900 

103 

73 

49 

56 

12 

7650 

130 

74 

56 

57 

12 

4501 

00 

75 

56 

58 

12 

2550 

130 

76 

58 

80 

8 

1560 

130 

77 

55 

59 

10 

3920 

130 

78 

59 

60 

10 

2820 

130 

74 

58 

61 

10 

4060 

130 

80 

60 

62 

8 

2960 

130 

81 

61 

63 

8 

2290 

130 

82 

33 

34 

8 

3330 

100 

83 

21 

70 

8 

440 

35 

84 

65 

75 

24 

70 

84 

85 

38 

74 

12 

1880 

102 

86 

29 

72 

16 

920 

69 

87 

41 

76 

18 

1550 

112 

88 

52 

78 

16 

962 

130 

84 

47 

79 

10 

540 

130 

40 

59 

80 

8 

1090 

130 

41 

62 

81 

8 

880 

130 

42 

63 

81 

8 

1000 

130 

100 

24 

90 

PRV  AT  85.0  ^si 

101 

16 

91 

PRV  AT  107.0  psi 
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Table  D2 


Node  Data  for  Fort  Hood  Main  Post 

Node 

Elevation 

Node  Demand 

Number 

(ft) 

(gpm) 

1 

860.0 

Reservoir 

2 

853.0 

0 

3 

853.0 

0 

4 

853.0 

0 

5 

853.0 

0 

6 

853.0 

0 

7 

853.0 

0 

8 

853.0 

0 

9 

853.0 

0 

10 

853.0 

0 

11 

853.0 

0 

12 

853.0 

0 

13 

853.0 

0 

15 

853.0 

0 

16 

854.0 

0 

17 

853.0 

0 

18 

859.0 

0 

19 

864.0 

29 

20 

906.0 

42 

21 

876.0 

16 

22 

877.0 

25 

23 

857.0 

31 

24 

877.0 

120 

25 

1049.0 

Tank 

26 

865.0 

0 

28 

902.0 

41 

29 

910.0 

73 

30 

911.0 

0 

31 

918.0 

100 

32 

1051.6 

Tank 

33 

892.0 

96 

34 

918.0 

2a 

36 

915.5 

90 

37 

940.0 

0 

38 

945.0 

0 

39 

1047.2 

Tank 

40 

953.0 

150 

4! 

936.0 

200 

42 

953.0 

110 

43 

911.0 

5a 

44 

950.0 

0 

45 

1051.3 

Tank 
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Table  D2  (Cont’d) 


Node 

Number 

Elevation 

(ft) 

Node  Demand 

(gpm) 

46 

912.0 

113 

47 

968.0 

106 

48 

955.0 

26 

49 

946.0 

111 

50 

1047.2 

Tank 

51 

952.0 

0 

52 

941.0 

88 

53 

941.0 

36 

54 

943.0 

291 

55 

898.0 

58 

56 

929.0 

117 

57 

1051.3 

Tank 

58 

909.0 

96 

59 

899.0 

54 

60 

847.0 

33 

61 

874.0 

64 

62 

870.0 

45 

63 

926.0 

52 

65 

877.0 

65 

70 

877.0 

0 

72 

906.0 

0 

74 

924.0 

0 

75 

912.0 

0 

76 

937.0 

0 

78 

930.0 

0 

79 

965.0 

0 

80 

888.0 

0 

81 

923.0 

0 

90 

877.0 

0 

91 

854.0 

0 
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USACERL  DISTRIBUTION 


Chief  of  Engineers 
ATTN:  CEHEC-IM-LH  (2) 

ATTN:  CEHECIM-LP  (2) 

ATTN:  CEMP-EB 
ATTN:  CECC-P 
ATTN:  CECW 
ATTN:  CECW-O 
ATTN:  CECW-P 
ATTN:  CECW-RR 
ATTN:  CEMP 
ATTN:  CEMP-C 
ATTN:  CEMP-E 
ATTN:  CERD 
ATr.N:  CERD-L 
ATTN:  CERD-C 
ATTN':  CERD-M 
ATTN:  CERM 
ATTN:  DAEN-ZCE 
ATTN:  DAEN-ZCt 
ATTN:  DAEN-ZCM 
ATTN:  DAEN-ZCZ 

CEUSC 

ATTN:  CEHSC-ZC  22060 
ATTN:  DET  IB  79906 
ATTN:  CEHSC-F  22060 
ATTN:  CEHSC-TT-F  22060 

US  Army  Engineer  Districts 
ATTN:  Library  (40) 

US  Army  Engr  Divisions 
ATTN:  Library  (14) 

US  Army  Europe 

ODCS/Engireer  09403 
ATTN:  AEAEN-FE 
ATTN:  AEAENODCS 
V  Corps 

ATTN:  DEH  (11) 

VII  Corps 

ATTN:  DEH  (16) 

21st  Support  Command 
ATTN:  DEH  (12) 

USA  Berlin 

ATTN:  DEH  (9) 

Allied  Command  Europe  (ACE) 

A1TN:  ACSGEB  09011 
ATTN:  SHIHB/Engineer  09055 
ATTN:  A  HUES  09168 
USASETAF 

ATTN':  AESE-EN-D  09019 
8th  USA,  Korea  (19) 

KOKA'S  Combined  Forces  Command  96301 
ATTN:  EUSA-HHC-CFC/Engr 

USA  Japan  (USARJ) 

ATTN:  DCSEN  96343 
ATTN:  Facilities  Engineer  96343 
ATTN:  DEH-Okinawt  96331 

Area  Engineer,  AEDC-Area  Office 
Arnold  Air  Force  Station,  TN  37389 


US  Military  Academy  10996 
ATTN:  Facilities  Engineer 
ATTN:  Dept  of  Geography  & 

Computer  Sciences 
ATTN:  MAEN-A 

AMC  -  Dir.,  Inst,,  &  Svcs. 

ATTN:  DEH  (22) 

DLA  ATTN:  DLA-W1  22304 

DNA  ATTN:  NADS  20305 

FORSCOM  (28) 

FORSCOM  Engineer,  ATTN:  Spt  Del.  15071 
ATTN:  DEH 

HSC 

Ft.  Sam  Houston  AMC  78234 
ATTN:  HSLO-F 
Filzsimons  AMC  80045 
ATTN:  HSHG-DEH 
Walter  Reed  AMC  20307 
ATTN:  Facilities  Engineer 

NARADCOM.  ATTN:  DRDNA-F  01760 

TARCOM,  Fac,  Div.  48090 

TRADOC  (19) 

HQ,  TRADOC.  ATTN:  ATEN-DEH  23651 
ATTN:  DEH 

TSARCOM.  ATTN:  STSAS-F  63120 
USA1S 

Fort  Huachuca  8561 3 

ATTN:  Facilities  Engineer  (3) 

Fort  Ritchie  21719 

WESTCOM 
Fort  Shaftcr  96858 
ATTN:  DEH 
ATTN:  APEN-A 

SHAPE  09055 

ATTN:  Survivability  Sect  CCB-OPS 
ATTN:  Infrastructure  Branch,  LANDA 

HQ  USEUCOM  09128 
ATTN:  ECJ  477-LOE 

Fort  Belvoir,  VA 

ATTN:  Australian  Liaison  Officer  22060 
ATTN:  Water  Resource  Center  22060 
ATTN:  Engr  Studies  Center  22060 
ATTN:  Engr  Topographic  Lab  22060 
ATTN:  ATZA-TE-SW  22060 
ATTN:  CECC-R  22060 

CECRL,  ATTN:  Library  03755 

HQ,  xvni  Airborne  Corps  and 
Ft.  Bragg  28307 

ATTN:  AFZA-DEH  EE 


NAVFAC 

ATTN:  Division  Offices  (1 1) 

ATTN:  Facilities  Engr  Cmd  (9) 

ATTN:  Naval  Public  Worits  Center  (9) 
ATTN:  Naval  Civil  Engr  Lab  (3) 

ATTN:  Naval  Constr  Battalion  Ctr  93043 

Engineering  Societies  Library 
New  York,  NY  10017 

National  Guard  Bureau  20310 
Installation  Division 

US  Government  Printing  Office  20401 
Receiving/Depository  Section  (2) 

US  Army  Env.  Hygiene  Agency 
ATTN:  HSHB-ME  21010 

Nat'l  Institute  of  Standards  &  Tech  20899 

Defense  Technical  Info.  Center  22304 
ATTN:  DTIC-FAB  (2) 
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416th  Engineer  Command  60623 
ATTN:  Facilities  Engineer 


Chanute  AFB,  IL  61868 
3345  CES/DE,  Stop  27 


